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ABSTRACT 

In recent years many theoretical and experimental 
examinations of the stability of chemically reacting systems 
have shown that these systems can display oscillatory behav- 
ior. In this work, a combined modelling and experimental 
investigation OlmmosScrinatory peneavioreinwtne. oxidation of CO 
over supported metal catalysts was performed. 

The initial modelling studies were of a general nature, 
and from them it was found that the use of simplifying 
assumptions in reactor models can result in transient and 
Steady-State violations of mass conservation. It was also 
found that complex periodic and chaotic behavior are possi- 
ble when a reactor model consists of more than two ordinary 
differential equations. 

From the experimental work, the detailed effects of 
operating conditions on oscillatory behavior in the oxida- 
Ciongon COMOVerpsupported Pt and iPts=Pd) catalysts were. deters 
mined. This experimental work enabled some insight to be 
Gained into the underlying mechanism responsible for oscil- 
latory behavior. 

Using the experimental insights, a specific model for 
the ene oxidation of CO was developed. For certain 
parameter values, this model successfully reproduced many of 
the fine structural details of the experimentally observed 


oscillations. 
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1. INTRODUCTION AND OBJECTIVES 


Ties LOcuct. On 
hiversced: Dyejenewcamby work vat Liljenroth [1.1] and van 

Heerden [71.2], a considerable amount of research has been 

directed at analyzing the stability of chemically reacting 

systems. Theoretical and experimental studies, summarized 
necevera® cvreview arcicies [1.3 - 91.7], have shown that 
under certain conditions, chemical and catalytic reaction 
systems can display oscillatory behavior. In a continuous 
flow reactor, oscillatory behavior is normally manifested 
through periodic changes in the reaction temperature and/or 
the concentrations of reacting species with respect to time. 

In heterogeneously catalyzed reactions oscillations in 

concentrations of adsorbed species and/or oscillations in 

the temperature of the catalyst's Surface can occur as well. 

These oscillations can occur without accompanying oscilla- 

ELONS ite the bulk fluid ‘phase, 

The recent interest in oscillatory phenomena in 
chemical reactors is due to several factors: 

1. The fundamental processes occurring in chemical reactors 
can be more fully understood through a knowledge of the 
mechanisms by which oscillatory behavior occurs. 

2. Large temperature fluctuations on catalyst surfaces can 
result in rapid catalyst deactivation. 


3. Large temperature fluctuations in a chemical reactor can 
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present safety problems, hence, a knowledge of how to 
predict, and suppress, this behavior is necessary for 
proper design and operation. 

4. Increased conversion of reactants to products, and 
isolation of low yield intermediates, can sometimes be 
obtained by operating.a chemical reactor in a cyclical 


Pashirons 


In the field of heterogeneously catalyzed reactions, 
reports of oscillatory behavior have been mainly restricted 
COmMOXxTGamronmreact10nS,— 6.07) Ox1dationm ot ammonia: (1.6.1; 
RyvowCCe rei .o 2 mis) seeCarbolr monox1de fal 2 7.1.91) and 
Vverocarbonshiie19e-9 1.21 le catalyzed by bulk (wire, gauze 
CEmLGit ) Oresupported platinum or nickel. Self-sustained 
ocillations have been observed in recycle, fixed-bed and 
single-pellet reactors under isothermal and nonisothermal 
conditions. Various causes for the oscillatory behavior 
have been postulated; these include competing mechanisms 
[1.13, 1.14], variations in activation energy with surface 
coverage [1.22, 1.23], non-reacting adsorbed species [1.24], 
presence of impurities [1.15], formation of surface oxides 
eae 1.25], variations in catalyst surface temperature 
UVecoliemand LOcadeneat) and mass poLanstermenteciuc melee 

Many mathematical models which display oscillatory 
behavior have been constructed using the above postulates. 
Unfortunately, a large chasm has separated the theoreticians 


from the experimentalists, as none of the models which have 
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appeared in the literature have been validated through 
Systematic experimentation. There are two main reasons for 
this: firstly, an experimental study which systematically 
examines the effects of the various operating conditions 
(temperature, concentration, space velocity, etc.) on oscil- 
latory behavior has yet to appear in the literature, and 
secondly, most models have been so reduced in complexity 
through the use of simplifying assumptions that only a 


tenuous connection with physical reality remains. 


1.2 Research Objectives 

As originally conceived, this research project was to 
encompass a general mathematical and experimental examina- 
tion of oscillatory behavior in heterogeneous catalytic 
reactors. Since a total description of this phenomena would 
indeed be a formidable task, a number of specific research 
areas, generally consistent with the overall objectives, 
were chosen. 

As a starting point, models which have appeared in the 
literature were examined with the particular aim of investi- 
gating the effects of simplifying assumptions, such as the 
equilibrium assumption. Also, in order to gain modelling 
experience, simple models displaying oscillatory behavior 
were developed. In parallel with this work, an experimental 
apparatus was constructed and oscillatory banavion Of carbon 


monoxide oxidation over supported metal catalysts was 
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examined. 

Building on the experience thus gained, a systematic 
experimental investigation of the effects of operating 
conditions on oscillatory behavior was performed. In 
parallel with the experimental study, a new model for carbon 
monoxide oxidation was developed. The new model was kept as 
physically realistic as possible, so that all parameters in 
the model had actual physical counterparts, and thus the 
model behavior could be directly compared with the actual 
experimental behavior. However, in a model, increased 
complexity is the price which is exacted for physical real- 
ism, and thus a total delimiting of the model behavior has 


not yet been achieved. 


1.3 Thesis Layout 

Since the research project was approached from several 
directions, this thesis is composed of several quite 
distinct though related pieces of work. For this reason, 
references are cited at the end of each chapter. 

In Chapter 2, a simple catalytic reaction iS examined 
SO as to develop the groundwork for determining the effects 
of simplifying assumptions on reaction models. Building on 
this development, the effects of the equilibrium assumption 
on a model of an oscillatory reaction are presented in 
Chapter 3, and in Chapter 4 the consequences of removing the 


equilibrium assumption, and going to a higher order system 
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(three ordinary differential equations), are examined. In 
Chapter 5, an interesting property of high order systems, 
namely chaos, 1S introduced, and is shown ne exist in a very 
Simple reaction model. 

The next two Chapters, 6 and 7, are devoted to a 
description of the experimental equipment, and a presenta- 
tion of experimental observations for carbon monoxide oxida- 
tion on supported metal catalysts. These experimental 
observations are then used in Chapter 8 to qualitatively 
examine the predictions of a new model for carbon monoxide 
oxidation. Finally in)Chapter 9 .possible' directions for 


future work are presented. 
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2. EFFECT OF SIMPLIFYING ASSUMPTIONS 


ON A NON-OSCILLATORY MODEL 


eee NcCroavuection 

In the modelling of heterogeneously catalyzed reac- 
tions, whether it be steady-state or dynamic modelling, it 
is common to begin with fairly complex models which consist 
of many equations. These large models result from the need 
to describe the many variables (bulk and surface phase 
temperatures, bulk and surface phase concentrations of the 
reactants and products, surface phase concentrations of 
intermediates, etc.) which are present in heterogeneous 
catalytic systems. 

Since the analytical or numerical solution of large 
systems of simultaneous non-linear algebraic or differential 
equations 1S a non-trivial task, it 1S common to reduce the 
Size of the system of equations by making simplifying 
assumptions. In the modelling of heterogeneously catalyzed 
reactions, the two most commonly used simplifying assump- 
tions are the steady-state assumption and the equilibrium 
assumption. 

These two assumptions have been used, usually without 
question, throughout virtually the entire history of hetero- 
geneous catalysis, as shown by examination of the early 
WOEKS [2 1,022.2 Nelms therfiedd. cist eusceonlyeasctor, Late: that 


any questioning of the validity of the assumptions [2.3], 
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and examination of their effects [2.4 - 2.6] has been 
performed. However, even in systems in which it can be 
shown that the surface reaction is the controlling step, the 
use of the equilibrium assumption in a model can have some 
unexpected effects. These insidious effects result from the 
very manner in which the equilibrium assumption is used to 
eliminate variables from a model. 

Rom nem@r Ol lOWInG yt wioll be wshownithat for a very 
Simple model, extremely different, and in some cases contra- 
dictory, results can be obtained depending on how the equi- 
librium assumption is used to eliminate variables from the 
model. For comparison purposes, the exact analytical solu- 
tion, and a solution using the steady-state hypothesis, are 


also presented. 


2.2 A Simple Model 
The irreversible isomerization of A will be used as an 
example. The overall reaction is 
A ——> B rare ie 


and the following reaction mechanism will be considered: 
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If the catalytic sites, X, are conserved,then mass balances 


ONAATaDXy anoexooiver 


CA xa Axe meet hee) AX aie eke [Ax G2 35) 
at 
Cpe eae ih eee kee DBKS 4 ke PAX (2°65) 
dt 
OU Ra teae Kee eee eA ak of Bal Ree kes PBX } Oe 287 ) 
dt 


where the normalized site conservation equation is given by: 
PARE EX ee exe =o 29S) 


In dimensionless form, Eqs. 2.5-2.7 are written as: 


opAeanks [Rh Sen etait eax) C2098 
ould 

Cnn) Bake Px aka, [BX aes Ax) (220) 
aT 

ke ee RO XA ce he AS oe Kee Perk SL Bx) 2m lal 
aT 

where, Kae ee eaal PR. (ear 

Ke, €2 "hah ke (22125) 

K, = Pes ks (PRC) 

Key ee heey ke (24 20) 

tT. =a (2242e) 


In the following it will be assumed that [A] and [B] in 


the gas phase are held constant. Conceptually, this could 
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be achieved by having only a very small amount of catalyst 
Liem Langemneactormmotmarternatively, it could be accom- 
plished by regulating the flows of A and B to and from the 


reactor. 


2.38 exacte solution 
Since Eq. 2.8 forces the sum of the derivatives to be 


G0UawEcOeZero, 1,65: 


GQpAx )£e d>EBXdt eal xd = 70 MEER 


aT dT aT 


then Eqs. 2.9-2.11 are not linearly independent, and any two 
of these equations are sufficient to completely describe the 
behavior of the surface concentrations. For example, if [xX] 
bowe Mimi nareca trom Eqs... 2.9 and 2.10 by using Bq..-2.8, then 


the following system of two ordinary differential equations 


results: 
CUA Sema wee ieK tKee) fe OR BX IT KG (Or 1a:) 

ar 
dl Bx] tae 1k oe Axaa—, (Ke FRE.) LBS] bes (2315) 

(ou 


From linear systems theory, the solution of Eqs. 2.14 and 
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aT lou 
PAS ook (Atk ee iaee™ aller tehe Ub+ hee). eee) 
a(a-b) bla-b) (2.16) 
aT Dr eulige aye 


+{(atK_,+K,)e [Conk oe thee } [AX], eee = BS, 
(a-b) (a-b) (a-b) 


wiscisi seoget tie 
Us es . « Wie -0'¥n{ at 99) Saat - 
| T 7 ag 4 7 - 

: 7 

et ag ag 
re Le m4 ‘ica , ton ate fh tee ghas Gare . y 
ue 

; 


fey og ly 4 .WJ¢ shh BO) Serie 6047 Se 
cs. Vet. Bal aia * pent > ed ro, fated 
hd ' ace . oot tet) Soderial fe ae 
wary | > neteye DnNetto sain 


: . . 7 See 1): i 7 
cars | Prat: 


o “Hye a“ 


Meas WEL I> U96cU.%& Bh), Sate coeza ye \eensy 


' Caee Ros 


aT lew 
BBX) Bemaik eke (ar iakoee (tiee  ) + {K,4K5 (b+14K.,)}(1-e ) 
ala-b bla-b 
at bT 
+{(at1+K,+K_,)e — (btitk,+K2,)e 1 BX],; 
a-b 
ak 9pt 
BiGioKe) hee) -ete)[ Ax}. (9517) 
a=b5 


where, 


ans (14K 4K oR, 4K.) [1 /TELR ORR RTT | 
2 (1+K,+K_,+K,+K_,)?2 
(2. ea) 


b = = UU a ol | TI Ea, 
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C2.upeD) 
Brom bos. 2. 16,ang 2.17, the values of [AX] and {BZ ] at a 


Steady-state are readily found to be: 


TAX] = =Ki(atKe,) + K,(btK..) 
ala-b) b(a-b) 
[AX] = eke (2.19) 
Kes at ae ema ( St Ree nh Ret Kone) 
[Bx] «= Rates Kier Koa) C220) 


KECIERS2) RBG. @) ake 


ThusyeEqs, «2al6—-2920 completely descruabe thetdynamic jand 


steady-state behavior of the reaction system. 
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2.4 Steady-State Assumption 

In the previous reaction system, the steady-state 
assumption could take the form of assuming that the concen- 
tration of either of surface species AX or BX is always at a 
steady-state. If, for example, surface species BX is always 
at a steady-state, then mathematically this is equivalent to 
SeutinogshG. 2. 10 equal Lorzerco, 

KACk) She UB eA X] (= 20 (So) 
DEMEG.)2.2i 2S USed 1m Conjunction with the site conserva- 
MuCnusequati1On) BEC necec, suhen tthe ysyscem of as; 2.9-2, i1.can 
be reduced to a single equation in one variable. 
Biiminatingetx’) and LBxX] from Ba. 9 gives: 


Gear hee ee PAK Ait Ke 2) +1 enh oe) CR +k. 3s C2: 67:38) 
aT K,+K_ > (KEK 
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eriaert a Rest te ee or ea) te Ke eee 
+[AX].e ee 
It should be noted that this result can also be 
obtained if Eqs. 2.8 and 2.21 are used to eliminate [AX] and 
[BX] from Eq. 2.11, as the same expression for [AX] always 
occurs irrespective of the solution path which is followed 
when the steady-state assumption is used. 
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Hien@ucamesteacvyestatembenavior, Also, 1£ KytK.5,>>1+K,+K.., 
then the dynamic behavicoreof Eq. 2.16 will be similar to 
CiatwonehGne2.23,easmenesconstantspa ana b in -thesexponen= 


tialmtenms of Eq: 2916" become: 


ane ~(K.4K..)[1- ATR OER CR RT| 
Z UK och Kase) = 
aie= = Kee eRe oe (tt Kee) (KotKk- 3} C2Ree:) 
hee hee 
Deeg has Rie etka Ith ot (at hee) AK OTEK WF Meee ea) 
K2ztK~, 

EcOmenos, 92.24 end 7.25 .1te1seseen ithat bis a much 

larger negative number than a when K,+K_,>>1+K,+K.,, thus 


bhesexponential terms containing b in Eq. 2.16-will rapidly 
approach zero, and hence the dynamic behavior of Eq. 2.16 
will be controlled by the exponential terms which contain a. 
This will result in the exponential terms of Eqs. 2.16 and 


Jeesebeing identical, hence the dynamic behavior of Eq’. 2.16 


Wore be asrmilaretoethatmorvngq.e 2523. 


2.5 Equilibrium Assumption 

The equilibrium assumption for this simple reaction 
system could take the form of assuming that either of AX or 
BX is at equilibrium with the gas phase. In order to main- 
tain consistency with the preceeding development for the 
steady-state hypothesis, it will be assumed that BX is 


always at equilibrium with the gas phase, which 
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mathematically is equivalent to: 

AP pe) a pe (2.26) 
BYBUSINGEECShe2e.26 ancy 2eo, akoasvpossible to) eliminate two 
mecvew lls aml hx | =eancCmuxImconcentrations from Bas. 
2.9-2.11, and thus reduce the system to a single equation in 
one variable. In the following it will be shown that 
extremely different solutions can be arrived at depending on 


which two variables are eliminated. 
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UeEIngeho.ececv LO eliminate [Xx ]) un ska.) 2.9 gives: 
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Comparing Eq. 2.29 with 2.23 reveals that the [Ax] from 


the equilibrium assumption will be similar to that from the 


steady-state assumption if K_,>>1. 
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2.5.2 Elimination of [xX] and [AX] from Eq. 2.10 
Eombinationfoferesy 2585and 2.26 cantbe used to yield 
bak @asma sounevLonson. | BX)", “such= that? 
CA eee em RE GK ths 20 7K SMES, 
Weing Equuziz4sto eliminate [xX], and Eq. 2.30 to eliminate 
RAS eC eel Overs: 


OR6s l= eee CB Ker Kee) (agua 
ag K, 


Solving hdeare.s 1 tom | Bx)-gives: 


au Wacko, ) Ke TRG Rae) 
BBX i= Ke ize ii PBX Iie (25 cee 
The behavior of [AX] can now be determined by substituting 
EGumceocmi tO a 2, 30 tOry1eld: 
pola Ket Keone Ke 
[AX] = [AxX].e estes, 

In order to remain consistent with the equilibrium assump- 
tion, Eq. 2.30 has also been used to convert [BX], to [AX], 
sg) Neus faces 

It is seen that the steady-state and dynamic behavior 
for [AX] predicted by Eq. 2.33 is extremely different from 
that given by EBq.e2.29, as Eq.) 2.33 predicts#a steady osuare 
valueloL LAX] .é6équaketo zero, which wlll only occum ingkg. 
Pood Sut Kiko, H0k However, even, 1D sthisecondicion, 1s met, 
the dynamic behavior of Eq. 2.33 will only be similar to 
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K_,.<<K, are met. It is interesting to note that these above 
conditions are necessary to secure dynamic and steady-state 
agreement between Eqs. 2.29 and 2.33 even though both equa- 
tions were derived from Eqs. 2.9-2.11 using the same equili- 


brium assumption. 


Zee i MinatlOneOLe Ax mands [Bx )) Prom Eq.. 2.5.11 
Supspucutionwol Ha 2 2601nto ha. 2.8 yields [AX) asa 
Punet. Oneo tas by 
PAs =a pe CRAKE 1') Kae C234) 
Using Eqs. 2.26 and 2.34 to respectively eliminate [BX] and 


UXmsin oho, e2.0isgrves: 
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aT Kea 
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The behavior of [AX] can now be determined by substituting 
EQwia. 36 ginto Eqti2.34 to gives 
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PRCOMDATUSCONMOrmEdS.. 2,37 and 2.29)shows that if K_,>>1 
and K_,>>K,, then similar steady-state and dynamic behavior 
Pore lA) awl ben preducted= py. bothsequations.-As-in-—the 
previous derivation, these conditions are necessary even 
ehoughebotn ot Eos.92.29 and 2,37 were derived from Eqs+ 


2.9-2.11 using the same equilibrium assumption. 


2.6 Conditions for Pairwise Equivalence 

LEstedeine Tables. 2.-kiand2ia2-are the wariouss conditions 
necessary for pairwise equivalent behavior between the 
expressions derived for [AX] using the exact solution, the 


steady-state approximation, and the equilibrium assumption. 
Table 2.1 Conditions for Pairwise Steady-State Equivalence 
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An inspection of the conditions listed in Table 2.2 


reveals some additional complications which result from the 


equilibrium assumption. Not only do the forms of Eqs. 2.29, 


Table 2.2 Conditions for Pairwise Dynamic Equivalence 
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2.7 Conclusions 

The preceeding derivations were not performed with the 
intention of illustrating that certain conditions are neces- 
Sary to secure agreement between an exact solution, and 
Solutions based on steady-state or equilibrium assumptions. 
This 1s intuitively obvious, and needs no proof. Rather, 
they were intended to illustrate that the steady-state 
approximation does not introduce any inconsistencies into 
the solution of a reaction model, while the equilibrium 
assumption does introduce inconsistencies. The inconsisten- 
cies introduced by the equilibrium assumption appear in the 
HOuM sO meciiterentepand contradictory, solutions to the reac— 
tion model depending on which variables are eliminated from 
the model by the use of this assumption. 

The question whose answer must now be attempted is: 
"What is the 'best' way of using the equilibrium assump- 
tion?". <A tentative answer to this question can be provided 
by examining the following equations which result from 


directly introducing the equilibrium assumption into Eqs. 


Dene Fania | We 
GHAS bh= th abs le Ke LAX) <— [Axa (22.38) 

aT 
Spex} =) Ax] (2.39) 

aT 
Gk =e ee! lek FAX] (200) 
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As was previously shown, the solution of Eq. 2.38 was 
Guyver OyebO. 2629 Manders gsathis equation, “of the solutions 
based on the equilibrium assumption, which most closely 
agrees with either the exact solution or the solution based 
on the steady-state approximation. A comparison of Eqs. 
Boon ce URWAEDRMEOS -299>2)iteshows thatthq. 2. 38hisithetonly 
eguation which has ee been explicitly altered by the equi- 
PIpGlumeessumptioge(anmimplicit alteration occurs in. the 
Site conservation equation). Hence, a very casual 
rule-of-thumb is to only use those equations which have been 
altered the 'smallest amount' by the application of the 
equilibrium assumption. 

BOrwtirs model, withthe analytical derivation 
Gomoreveq asic cam be contidentiy said that the implicit 
alteration of one term in Eq. 2.38 is a smaller alteration 
than the total elimination of two terms which has occurred 
in each of Eqs. 2.39 and 2.40. However, in general there is 
no way of determining if one type of alteration is smaller 


than another. 


2.8 Nomenclature 


Roman 
a,b Gonstantsvinethe analytical solutionsol Bos. 2.142.515 
[A] gas phase concentration of A 


[AX] normalized surface phase concentration of A 
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gas phase concentration of B 

normalized surface concentration of B 

rate constant for the i-th forward reaction 
Parenconstantemoreene Iath forward reaction 
dimensionless rate constants for the forward reactions 
dimensionless rate constants for the forward reactions 
time 

dimensionless time 


normalized active Site surface concentration 
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3. EFFECT OF THE EQUILIBRIUM ASSUMPTION ON A 


MODEL FOR AN OSCILLATORY REACTION 


Sra en troduct ton 

In the previous chapter it was seen that the use of the 
equilibrium assumption can have unexpected effects on both 
the steady-state and the dynamic behavior of a reaction 
model. Since these effects were present in even a very 
Simple model, they should also be present in the more 
complex models which are used to describe oscillatory 
behavior in chemical reactions. In order to determine the 
effects of the equilibrium assumption on a model of an 
oscillatory reaction it was decided to examine a model 


presented by Eigenberger [3.1]. 


3.2 Eigenberger's Model 

Eigenberger proposed that oscillatory behavior in cata- 
lytic reactions could be caused by the existence of an inert 
Surface species which alternately acts as a source of, anda 
Sinks £Orecatalytic sites, For the overall, reaction 


2A + B —» 2C, Eigenberger proposed the following mechanism: 


B + 2X BXX (Fo 2) 
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OARS kee se ioe 4X 33) 


Berax BX (3e4s) 


The adsorbed species BX is unreactive in this mechanism and 
actS aS a capacitance which makes oscillations possible. 

In several papers [3.1 - 3.3] devoted to describing 
this model, it was always assumed that the BXX species was 
at equilibrium with the gas phase. In this way, a reaction 
model describing this system can be reduced to a system of 
two ordinary differential equations, as two of [AX], [BX], 
[BXX] and [X] can be eliminated through the use of the equi- 
librium assumption, and the site conservation equation. 
However, based on the work in the previous chapter, the 
question must arise as to which two variables should be 
eliminated. 

A general model for this system based on mass balances 


on the surface species 1S as given below: 


CAS Ie aka AX Be Pe A Be) 3 5 

at 

aLBx =k EB lo - k  bB] (3.6) 

Nokes 

GPBXX |. = kh, UB bik = ok 2 DBRS ek PA See (Sei) 
at 


Gheee= kh RAR eae VA) LS) ce ee Ree Ree Ck) 
dt 


tol ehe ke fer [ BULK IES 4h hak Bxx] (3.8) 
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where it is assumed that [A] and [B] in the gas phase are 
held constant by some means or other, and where the norma- 
lized site conservation equation is given by: 

Pex cee x tech ea, (xe 8H (39) 


In dimensionless form, this system of equations is given by: 


CA =e. | Soe ees Aero AX Pe Bex] ¢G3..10) 
oul 

CUBA Patel Xl E xX) (3.14) 
aT 

GU k Rea ho eee aK ee | BEX] = OK PAX TBR] (Soe) 
aT 
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K, = k,L[B]/k_, (3.4c) 
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and by using the normalized site conservation equation, 

[AX] pet BX ere a= 7 CRG 
In the site conservation equation, Eigenberger assumed that 
[BXX] was always very small, and thus was insignificant 
compared to the other terms. However, the development in 
the preceeding chapter has shown that the model performance 
is affected by the choice of variables which are eliminated 
through the use of the simplifying assumptions. Thus, it is 
necessary to decide which two. of) [Ax], [Bx], [x] or [Bxx] 
should be eliminated through the use of Eqs. 3.15 and 3.16. 
This decision can be partially resolved by examining the 
resultant torms Of Eqs. 3.10-3.13 if the equilibrium assump- 


tion 1s assumed to hold, i.e. 


GRACO N= ERO een FRA Se eek, (AX } be xx] C4 
aT 

AUB xa he LX ee LEX (2a te) 
aT 

eek Wee hk BX S| (oO) 
aT 


dik) =sKe, pAR ye Rex) + (BX) -Sk pk re eke AR Te eo) 


Clearly, Eq. 3.19 should not be used; as use of this vequa— 
tion will have one of two results; either [BXX] will be 
driven to a stable steady-state at zero, or some of the 
other concentrations will be driven into physically impos- 
Sible regions, i.e. outside of the interval (0,1). Thus, as 


was done by Eigenberger, Eq. 3.15 will be used to eliminate 


oe eae nk 
ry sdvienta , ah ania ig RTS aisle witd ni 
hf ~ : “! Bey AD SAA ell hee | [ew 424 a diene (owas 


aes j ; 
me | i. 9 rears Ve 
od Qi eseeonag, 


7 


: Sn Bi @hw we , > bia eu ‘ee 


‘ 


1 ps aides il jy 


og - 


Zu 


[BXX], with the following set of equations resulting: 


Cpe aaaae he | Me Ree PAS | ee Ko TAX |e xX) 2 (3824) 
aT 

CPB eek xe orl (2505 
aT 


CC ae ee Aree ee eR Xela Kp X ee CK TAN Re Ex Ie (3 23) 
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where, Keer = 21 hee (324) 
However, the choice is still unclear as to which of the 
remaining three variables should be eliminated by Eq. 3.16. 
In the following, the three possible different systems, each 
consisting of two ordinary differential equations, will be 
examined. In particular, the steady-state and dynamic 


differences between these systems will be illustrated. 


3.3 Elimination of [X] — System #1 
Pe 2G. 3.16 2s8used!/toseliminate [xX] efromi Bas: 
S.2t-2a2a, then thelfoblowingrsystem, to be known as system 


Pi eeresuLlts.: 


GUAM S=Ah AG LAM a Ox yee eh AM) 
aT 
= Kee LAX) ?( 4—[Ax |=[ Bx (3.25) 
abba =AKmC ipa oro xX) vee oa! (3.26) 
aT 


In order to better understand the behavior of this system, 


both a Uniqueness, and a stability analysis will be 
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performed. 


3.3.1 Uniqueness Analysis 
A uniqueness analysis consists of determining when the 

Ssystemmiofekhqss 3225 and 3.226 has*a unique ,steady-—state.This 

can be determined by performing the following steps: 

1. The left-hand sides of the differential equations are 
set equal to zero. 

2. The resulting algebraic equations are combined to obtain 
a Single algebraic equation in one variable. 

Se) Therconditions forawhich#this algebraic equation has 7a 
unique zero are determined. 

Setting the left-hand sides of Eqs. 3.25 and 3.26 equal to 


zero gives: 


Kee (PASI EBX) Rev bAs 12 = Roe = (AR FH Bey) He 20 (ree ys 


Veg ee oP hea ee (3428) 


Souvinomtomss 2onrores| BXJip) and Subsceltucing thisminto sq. 


3.27 gives upon rearrangement: 
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Applying Descartes Rule of “Signs towEq. 3.29 reveals that it 
has one or three positive real roots, and one negative real 
root. Since only roots in the interval Os([AxX]<? are physi- 
cally reasonable, it is necessary to check the number of 


positive roots which are outside of this range. This is 
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donesbyesubstituting [AXJ=Zt1einto EBq.93729 tonget: 

ee ea om er eh lth) oe Kea ith) =) 0 03.30) 
ee K>3 

An application of Descartes Rule of Signs to Eq. 3.30 shows 

that there are no roots in the region Z>0, hence there are 

no roots in the region [AX]>1. Thus, the region 0<[AX]<1 

can contain either one or three positive real roots. 

The conditions under which only a single positive real 
root will exist can be determined by examining the discri- 
Mitianteor Eq. 3.29. “For the general quartic ‘equation 

ue ige Japa a8 Mebes tys tebe tim Kova me (sta) 


the discriminant is given by [3.4] 
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where, fa: Sr derek, SS ee es (Sea c) 
and Ce eee DCO cme EOCe 25 mec ec 2g (3530) 


The type of roots of the quartic equation can be determined 


from the sign of the discriminant as shown in Table 3.1 


Table 3.1 Sign of the Discriminant and Quartic Root Types 
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Single poSitive real root is that the discriminant of Eq. 
3.29 should always be negative. Of greater interest is the 
Situation where the discriminant is equal to zero, as this 
delimits the boundary between the regions of unique 
steady-states and multiple steady-states. 

Phesc1scrimingne -ot Eq. 3,295 1s" found by substituting 
EnemCOGEM@ICRentS Of ong ms. 29 aro ko? 3.21 to give after 


considerable manipulation: 


hea 27 OP +O) + 4 (PO) (P2+34PO+0" ) — 16P0O (3.32a) 
where, Dare ale K oo) aK. 3.325) 
Ome Koala thei ceek 2 & (27220) 


Setting the discriminant equal to zero and rearranging the 


terms gives: 
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Thus, Eq. 3.33 can be used to define the boundary between 
the regions of uniqueness and multiplicity. The easiest way 
tomuse Eq, 3.33 018 teosmnotice that tt 1s 7a qQuadtetze an terme 
of Kip) if Ka haa and K. are specified, then the values 
(if any),ot Kz; fer the uniqueness > multiplici ct yaboundauy 
can be calculated. If the solut#on of Eq. 3.33 results in 
complex conjugates for K,, then uniqueness 1S guaranteed. 
Shown in Figure 3.1 are typical regions of uniqueness and 
multiplicity in the K,, - K, parameter space, where K,=48 


and K.,=2 were hela constant. 
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Figure 3.1 System #1 - Regions of Uniqueness and Multipli- 
CitveaSvasrunction of) Ks tandeky) (K,=48, K_\=2) 
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go3.e StabidityPAnalysis 

A stability analysis for the system of Eqs. 3.25 and 
3.26 can be performed by examining the signs of the trace 
and the determinant of the Jacobian matrix, where for this 


System the Jacobian matrix iS given by: 


7 9 ayer rigs eA CR iA kl le PAR Ba) 
-2K,;,[AX](1-[AX]-[BxX])? -K, 
TAL Ax A iL AMI at Bx) 

pee (3534) 
=k at hun Be 


A steady-state will be asymptotically stable if all of 
the eigenvalues (roots of the characteristic equation) have 
negative real parts. For a system of two ordinary differen- 
tial equations, the two eigenvalues will have negative real 
parts if the trace of the Jacobian is negative, and the 
determinant positive, where the trace and the determinant 
are evaluated at the steady-state. If either of these two 


conditions are violated, then the steady-state is unstable. 
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From the Jacobian matrix the trace 1S given by: 
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The number of unknowns in Eq. 3.35 can be reduced by two 
through the use of Eqs. 3.27 and 3.28. These equations are 
most conveniently used to eliminate [BX] and K.;, as solving 
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PBR) = KCI VAx 107 C1+K, ) (3.36) 
Welch mUDOMESUDS CI CULwOne Into mG. Se2/, ana solving for KZ; 
yields: 


Kavos ik kee LAR oho AS lth, ) (93. 3a) 
PN Wa fg) ae 


PiSeGEangenos, 3.36 >and 3.137 AMmtethot. 3.35 Gives: 


ee nee ahem Ax ie 2Re tte R TARY (937.38) 
(1+K,) [AX] 1-[AX] 
In order to determine the boundary at which the trace 
changes sign, Eq. 3.38 is set equal to zero and rearranged 


to give: 
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EO mma ese meOtevatlestifOluh,, hoy and. Kay BOs 3.08 Can be 
solved for two values of [AX] at which the trace is equal to 
zero. These values can then be used in Eq. 3.37 to deter- 
mine the corresponding values of K., for which the trace is 
euuaEcOezeron me Utrtnen values of VAX |yirom Eq 43. 39eare 
complex conjugates, then the trace is always negative, and 
for the selected values of K,, K.; and K, it cannot change 
Sichecor, any Kos. shown in Pigureémo. 2eare the pointssinehae 
- K, parameter space at which the trace of the Jacobian is 
edualeto Zero, ewhere the Values Obs manda kyaware these 


which were previously used in constructing Figure 3.1. 
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Figure 3.2 System #1 - Trace Condition Stability Regions as 
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Boe cn ee Determinant Condition 

From Eq. 3.34 the determinant of the Jacobian is given 
by: 

Dem ome Ga ke Kem k | PAXM 1 (ARN - (Bx dss 


+ Kia ORG [AX 18 Ot=-([Ax] = BK) (3540) 
SUbSUECUE IGehdS 6S. 36uald 44/8 intonned. 3.40 qives: 


Det omaha mck GltKe Anes RE, (1th, ) = 3K, (341) 


[AX] (= TAxd 


Setting Eq. 3.41 equal to zero and rearranging gives: 
he eee Ka) Ak wee OKO Ke ie) PAR] Son Seo (3.42) 


BOmanve Given set ofht@parametens Ky, Kl, and Ka, Bq.3.42 can 
be used to determine two values of [AX] for which the deter- 
minant of the Jacobian is equal to zero. These values can 
then be used in Eq. 3.37 to determine the corresponding 
Values soto K5, for whach the determinant is equal to zero. 
If complex conjugates result from the use of Eq. 3.42, then 
for the selected values of K,, K., and K,, the determinant 
LSepost@uivel forvablevalues of Kz.. 

When the above procedure was numerically implemented, 
Me wase foundsthat the locus of points at which thesdecer. 
minant was equal to zero, was identical to the locus of 
points which separate the regions of uniqueness and multi- 
plicity. Hence, the curves in Figure 3.1 represent the 
locus of points for which the determinant is equal to zero, 


as well as separating the regions of uniqueness and 
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Mb plicity. eThiceresult, 1S consistent with the 
theoretical findings of Gavalas [3.5], who showed that for 
systems obeying conservation equations, the number of 
Steady-states will be equal to 2mt+1 (m=0,1,2,etc.) where at 
least m of the steady-states will be unstable. 

In order to show the overall multiplicity and stability 
properties of System #1, a combination of Figures 3.1 and 


3.2 as shown in Figure 3.3 is necessary. 


3.4 Elimination of [AX] — System #2 

If Eq. 3.16 is used to eliminate [AX] from Eqs. 
3.21-3.23, then the following system, to be known as System 
#2, -results: 
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aT 


+ 2K,,(1-[BX]-[X])?[x]? (2a) 


hE x fei Kable BX] (3.42) 
aT 


3.4.1 Uniqueness Analysis 
Setting the left-hand sides of Eqs. 3.43 and 3.44 equal 


to zero gives: 
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Figure 3.3 System #1 - Combined Stability and Uniqueness 
Regions as a Function of K,,; and K, 
(K,=48, K_,=2) 
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SOs Cen mo NT OMGOcMEB Srsanc SUDSTItUtIng into Eq. 3.45 


gives upon rearrangement: 


Ok ite ket itKy Tk) * os 2Kh5 ba? 


=e Ket K oa Veh Pk) te Kiyo 0 (3047) 


However, at a steady-state, combination of Eqs. 3.16 and 
3.43 yields a relationship between [X] and [AX] such that: 
le Air iARI 7 hei thie) (3.48) 


Substitution of Eq. 3.48 into 3.47 gives upon rearrangement: 


PA he 2 pak) eee LAX 2 (3.49) 


A Goa) PCL eRe ee te Ke Cth Ef 2K>3 = 0 


A comparison of Eqs. 3.49 and 3.29 shows that these 
equations are virtually identical, with the only difference 
being that the—Ky. terms in Eq. 3.29 are multiplied. by a 
Pectorecm .wo- tnerores. 49eee Thus) Bthesresults toteagunigue— 
Becomdialysas tor System #2 will be rdentrcal to that) whieh 
was performed for System #1, with the only change being that 
the values of K,, necessary for multiplicity will be 
one-half of those previously determined for System #1. If 
the values on the K,,; axiS in Figure 3.1 are divided by a 
factor of two, then the modified figure would describe the 


uniqueness - multiplicity regions for System #2. 
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3.4.2 Stability Analysis 
From Eqs. 3.43 and 3.44, the Jacobian matrix for System 


wee ce given by: 


Aker PB ll) ol ol aR ake, he a5 NG 

Sheek ian Bore a ake poh erael NB Se) Ls bal ole 
J= e250 
Kee i 


If the Jacobian is written in terms of [AX] and [BX] through 


the use of Eqs. 3.48 and 3.16, then it becomes: 


Ae AR he Cie PAR Te RB XL) ae ee, 
-4K,,[AX](i-[AX]-[Bx])? -~4K,,(AX](1-[AX]-[BX])? 
J= Rs Ke ai Cao 41)) 
Ke oe | 


344.2.) Sl race t€ond tron 


Fromelq..3 15 ithe ‘trace sofathewJacobtanvis given sby: 


Pee eee he a a Ke AX uta PAX |B s 4d) 


ee A Rel Cle bA | ALBAN) (3252) 


Compare nguha. so. S2ewith 93.85 Shows! that the trace cta the 
Jacobian Lor Systemuf2 1S identical withethat cle Systemasu, 
except that all the K,, terms have been multiplied by a 
factor of two. ithus, if the values vonmestnewhs  axticsos 
Figure 3.2 are reduced by a factor of two, then this figure 


can also be used for determining the stability of System #2. 
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Ste ce Dererminant Condition 

From Eq. 3.51 the determinant of the Jacobian is given 
by: 

Dee) am 4 Ke) Reena ke AXA PAX] -f BRI) 23 


See UA Sl au AKA eal) (2m5o) 


Comparing Eq. 3.53 with 3.40 shows that the determinant of 
Phemeaccoren wor System #2) 1s identical with that of System 
#1, except that all of the K,, terms have been multiplied by 
SmEoOCtCOmnOL EwoOw sence wif Figure 3. 3 l1s modi fied «by 
dividingscne values on the K,, axis by a factor of two; then 
it can be used to depict those points at which the deter- 


minant is equal to zero for System #2. 


3.5 Elimination of [BX] — System #3 
MoniGmed eo aloe lose Iminate a BX) Morom EQS a Shei ss 


Gives System #3: 


Clee ere ee ee Ai he hax] 2 [x] (3554) 
oly 
Ce ee eer aA, ee a od rt nd a Os |e me OA] 
ar 
Tee? Ree eel ee (3.55) 


3.5.1 Uniqueness Analysis 
Settinglithe leftahandisideskot (hast eon odiand 2455eequal 


toezerocgives: 


Roto T ane PAX am Ki Aoi ne ee ae (3.56) 
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Beh |e kes PAR x et dO< P(X) Soraxy 


che UAK lek | 2yewe( ee 5H) 


Combai nom cews >omandmwoy7eGtves [xX |einy terms, of [AX] “such 
that 

xs eee eA ether eae Chet Kom 1.) (3556) 
SHDSeL Cur TOnmOL EG. oeoG incom woo yields wpon rearrange- 


ments 


Kereta PA emcee i tke se) PAX | eer Ko5 CAR.) 2 


—~ {K,+K_,(14+K,)}(K,-K,.-1) [AX] + K,(K,-K,-1) = 0 (3.59) 


Theougn che use of Descartes Rule of Signs, it is seen that 
Eq. 3.59 can possibly have zero, two or four positive real 
BOObSmi fuk, =Kk,tt,..On ab can have one or three positive real 
roots if K,<K,+i. As previously explained, it is necessary 
to determine how many of the positive roots are in the 

Bree hVa cS PAX <u sSUDSLT tubing VAR J=Z+h into Eq 4 3.59 


gives: 


peer te ots et ee ee ae Kit) a Ke eet eG Ke a eee 
3 Te ee ere 2 esi) ad Reet Roe ae er hay a ya 


PUA Kee Kh oer (al Fhe, (Kae eh ea) (3.60) 


Lieke<K ets Chenwalde Or ethe, cocii i cvent sane, as .oUma Ge 
positive, and hence there will be no roots in the region 2>0 
GRAX)> h) SheThusyethetreqionu0s[ Ax) <tewiteconcainmonesor 
three positive real roots when K,<K,+1. However, if 


K,>K,+1, then depending on the value of K., there can 
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POsscr bly beezeropmonenor tworrootsein thesregion [Ax]>1. 

Obviously, these results indicate that something is 
drastically wrong with System #3, as roots in the region 
[AX]>1 are physically impossible. Also, as stated 
previously, Gavalas has shown that systems which obey 
conservation principles must have an odd number of 
Steady-states. Since, for certain conditions, System #3 can 
have an even number of steady-states, then the conclusion 
which must be reached is that the use of the simplifying 
assumptions has resulted in a violation of mass conserva- 
ELON. 

In order to further illustrate some of the unusual 
features of System #3, the boundaries separating the regions 
of differing numbers of steady-states are determined by 
examining therdiseriminant of Eq. 3.59. Substituting the 
SQCPELCPENUSPOrRrhGs Geoldeintosko. 3<3h“yieldspafter consid- 


erable manipulation: 


eee Give hie Phen Onee ae 2 LR PA Oe 


PCS Uh an oO. Re OO (29GK7 Ry OF TOCK oO ©) PAKe a ko. Gua) 


where, Pa=mK an Knee | (3.61b) 
Wey oe ye co 3.6 103) 
Re=ekoetekhan(a tke) (361d) 


Setting the discriminant equal to zero and rearranging the 


terms gives: 
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Hotere ee eee oa Ke Re Omar W96K27RO7 "+ 62K202)K.% 


Toe PORs = 0 (ES) 


SunNCemud ss O2mtsSwaecuagCratic equation in iLerms-of K% a 
PueCHmllIh we ih Seana k, are specitied, the valuessotwks., for 
which the discriminant is equal to zero are readily deter- 
mined. These boundaries are shown in Figure 3.4 which has 
been constructed using the same values of K, and K_, as in 
Picurems |. lt Seinteresting, to note) that in Figure 3,4 
there 1S a region of parameter space for which System #3 
does not have any real steady-states. This is obviously 


physically impossible. 


3.6 Comparison of Systems #1, #2, and #3 

MASOUGne cnhewuseeoltehiroures 37.5 and (3.4 1 “1s possible 
POuchoose parameter yalues for K,4 and@k, such that inter- 
esting differences between the three systems are illus- 
trated. In all cases, the parameters K,=48 and K._,=2 are 
used, and the computer programs used for the integrations 
are included in Appendix A. 

For the parameter values K,;=20000 and K,=50 an 
examination of Figures 3.3 and 3.4 shows that all three 
systems possess a Single, stable, steady-state. However, as 
shown in Figure 3.5, the location of the steady-states 
@iffer significantly from each other. Thus, despite the 
three systems being qualitatively similar, they are quanti- 
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Figure 3.4 System #3 - Regions of Differing Numbers of 
Steady-States (K,=48, K_,=2) 
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[BX] [AX] 


Fractional Surface Coverages 


[X] 


Dimensionless Time, T 


Figure 3.5 Steady-State Differences Between Systems #1, #2, 
and #3 (K,=48, K.,=2, K,;=20000, K,=50) 
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behavior is analogous to that depicted in the previous 
chapter, where it was analytically shown that the location 
of a steady-state could change depending on the manner in 
which the equilibrium assumption was used to simplify a 
model. 

Perhaps a more interesting situation is shown in Figure 
3.6 where the parameters K,,;=1250 and K,=5 were used with 
the previously stated values of K, and K_,. From the 
uniqueness and stability analysis it is found that System #1 
has a single, stable, steady-state; System #2 has a single, 
unstable, steady-state; and System #3 does not have any real 
Steady-States. From Figure 3.6 it is seen that these 
differences result in dramatic qualitative, as well as quan- 
titative, differnces between the three systems, as System #1 
is stable, System #2 oscillates, and System #3 'explodes'. 

It 1S the behavior of System #3 which is particularly 
interesting, as if the integration 1S continued past that 
shown in Figure 3.6, then the value of [X] approaches posi- 
tive infinity, and that of [BX] approaches negative 
intaneey. loves behaviorliys. a directs vigolatiGn obemass 
conservation as all the concentrations have been normalized 
torbe in the range of zero to one. 

While it is apparent that these effects are due to the 
use of simplifying assumptions in the derivation of the 
three systems, it must be remembered that two simplifying 
assumptions were used. These assumptions were the 


assumption of equilibrium given by Eq. 3.15, 
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KX ieawk eb BAX | ae, 
and the assumption of a negligible concentration of [Bxx], 
which affected the site conservation equation (Eq. 3.9) to 
GUV Gs EC satsi8) OF 

PAX PA ePB Rete x li | (3.406) 
The previously shown violations of mass conservation could 
be due to the effects of either, or both, of these Jars vitae 
tions. In the following section, the [BXX] term will not be 
eliminated from the site conservation equation, and the 


resulting effect of this on the three systems will be shown. 


3.7 Isolation of the Equilibrium Assumption 
If [BXX] is not assumed to be negligible, then the 
normalized site conservation equation is given by: 

PA OX jee eet ODE Xx ee 03.363 ) 
However, if the BXX species 1S assumed to be at equilibrium 
with the gas phase, then [BXX] in Eq. 3.63 can be eliminated 
ChrOuUCGMmeEenemUSe aCletbg aus. HS Eto hyie Hd: 

PAS At EEX Weee (oa), RAK kis / Kat. 37.64) 
Equation 3.64 can now be used to eliminate one of [AX], [BX] 
or [X] from Eqs. 3.21-3.23 to arrive at three new systems, 


each consisting of two ordinary differential equations. 
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3.7.1 Elimination of [AX] — System #1A 
KyOnelaut >. ot metne pnystically realistic (0s[X)<10) solu- 


BeOne Cle Sl alSyGlvenDy: 


[x] = Ka [1+ / + 8K,(1-[AX]-[BxX (3.65) 
4k, he 


thts pe oyS een 7 JA iSugquven by pds. 3.21 and 3.22, where [xX] 
1s? determined, from Eq. 3.65. System #1A is virtually ‘ident- 
ical to System 41, with .the-exception--that only~onexsaimpli- 
fying assumption, the equilibrium assumption, is used in 
System #1A. The removal of the assumption that [BXxX] is 
negligible results in System #1A containing an additional 
Harameter Group, KoyKue, which: does not explicitly occur in 
System #1. In comparing the behavior of Systems #1 and #1A, 
it is necessary to specify a value of K,/K_,. The value of 
K,/K., should be selected such that [BXX] is very small in 
System #1A, thus the effects of the equilibrium assumption 
can be isolated. 

Shown in Figure 3.7 are the results of integrations of 
System #1A for several values of K,/K_.,, where the other 
parameters have been held constant at the values used in 
Figure a6. Whe initial values: (AX) .=0. 2) e1Bx),.-00c, 
[X],=0, and [BXX],=0 have been used for each integration. 
These values were chosen to maintain consistency with Figure 
3.6. Lteshould be noted that) £ a) x) 40 gic isedes then 
[BXX],#0 is necessary for the initial conditions to be 


consistent with the equilibrium assumption. 
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Fractional Surface Coverages 
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Figure a System #1A 
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Peon es GUY Cusamaicemseseen that tf Ko/Kio5 016 small; 
then the values of [AX] (also [BX] and [X] which are not 
shown) from integrations of System #1A are essentially 
Mleimrcalm Omeliatmercolmoyvercemnaish,/K.3 = 0)..- This, 15 tro 
be expected, as through the use of Eq. 3.16 it is seen that 
See valves  OLluKe A wi lerorcestBxx |] to be small relative 
to [X]. Hence, the assumption used in System #1, that [BXx] 
1s negligible relative to the other concentrations, will be 
Sauistied when K>/K_, 1S Small, and under these conditions 


System #1 and #1A are identical. 


3.7.2 Elimination of [AX] — System #2A 

From Bq. 3.64, [AX] can be written in terms of [BX] and 
Lxeletewoi ye 

(oe) ieee me ge Slee A ee Xr) a Kins (3. 66) 
Shes, system g@ZA 1S given by Eqs. 3.22 and 3.23, where? [Ax] 
PSmCeterminederrom Eq. 3.66. Shown in Figure 3.8 sare the 
results of integrations of System #2A for several values of 
K,/K.2, where the other parameters have been held constant 
at the values used in Figure 3.6, and the initial conditions 
were the same as those chosen for Figure 3.7. It is seen 
that when K,/K_, 1S small, then the behavior of System #2A 
VSavirtuallyuidentical cto that. Ohvoystenesc. iC mie We meas 
increased from that used in Figure 328(b)y sthen the unstable 
Steady-state becomes a stable steady-state as seen in Figure 
3.0(C) “This ocectrs due to the sSurftace coverage of the BXxX 


species becoming significant at large values of K.,/K-.. 
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However as Shownin migurne 3.610), 1f Ki/Ke; is further 
increased some additional interesting, but unexpected, 
behavior occurs, 

In Figure 3.8(d) it is seen that System #2A can display 
a dynamic violation of mass conservation when K,/K_, is 
Vamoe;, sas the jconcentration of AX decreases from its initial 
value of 0.20 and passes through negative values before 
finally reaching its steady-state. Thus, steady-state 
conformity to the principle of mass conservation does not 
necessarily imply the same for the system dynamical behav- 


Or. 


3.7.3 Elimination of [BX] — System #3A 

Imma Manner Simpler to the above, [BX] in-terms of -[Ax] 
and [X] is given by: 

eer ee ean Ace eC eee Re a) 2) Kees tse) 
Thue wSyStemasen, iS given by Eqs. 3.21 and 3.23, where |[ Bx] 
Pseoetermined wrom Eq. 3.67. Shown in Figure 3.9 are ghe 
results of integrations of System #3A for several values of 
K,/K.2, where the other parameters and the initial condi- 
Bions are theysamelas Those uscd mer oures Osea cape l 
is seen that removing the assumption that [BXX] is neglig- 
ible has not improved the performance of System #3A relative 
to System #3. 

For small values of K,/K., the behavior of System #3A 
TSRL Cent Ca lelCOmChaLMOc eSYSLEM #3. IihBthdmeonemviOlLaclons. Of 


mass conservation which occur in System #3 are also found to 
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Sccun unesystemsA3Awm Sincesonlysthesequiliibraum assumption 
has been used in deriving System #3A, then it can be 
concluded that these violations of mass conservation are a 


direct result of the use of the equilibrium assumption. 


3, 8Conclusions 

The use of the equilibrium assumption in a model of an 
oscillatory reaction has resulted in three simplified 
models. Despite the intuitive belief that the predictions 
of these models should be identical, it has been seen that 
Quantitative and qualitative differences exist between the 
three models. These differences occur in both the steady- 
State and the dynamic behavior of the models. In particu- 
lar, the use of the equilibrium assumption can result in 


dynamic and steady-state violations of mass conservation. 


3.9 Nomenclature 


[A] (gas phase concentration of A 

[AX] normalized surface phase concentration of A 

[B} gas phase concentration of B 

[BX] normalized surface concentration of reactive B 
[BXX] normalized surface concentration of unreactive B 
b constant in the general quartic equation 


e constant in the general quartic equation 
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a constant in the general quartic equation 

e constant in the general quartic equation 

J Jacobian matrix 

ee rate constant for the i-th forward reaction 

Kia EacehconSeant.slomecnem@arth;rorwardireact pon 

K; dimensionless rate constants for the forward reactions 
Ie dimensionless rate constants for the forward reactions 
Kos combined equilibrium and reaction rate constant 


P,p grouping of constants used in the discriminant 


OF grouping of constants used in the discriminant 


R grouping of constants used in the discriminant 

t time 

a6 dimensionless time 

Variable in the general quartic equation 

[x] normalized active site surface concentration 

Z origin shift variable used in the uniqueness analysis 
Greek 

A discriminant for the general quadratic equation 
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4, EFFECT OF THE REMOVAL OF THE EQUILIBRIUM ASSUMPTION 


4.1 Introduction 

In the preceeding chapter, it was seen that the use of 
the equilibrium assumption in a dynamic model could have 
unexpected, and unwanted, effects. In particular, it was 
shown that the manner in which the equilibrium assumption 
was used could have a dramatic effect on the steady-state 
and dynamic behavior of the resulting model. These effects 
could in some situations be manifested through violations of 
mass conservation. In this chapter the further effects due 


to the removal of the equilibrium assumption are examined. 


4.2 General Model 
In the previous chapter it was shown that the use of 
the mechanism proposed by Eigenberger resulted in a model 


composed of the following equations: 


OC CAxs| R= (Fenix hehe Ae ere AR See] (4.43 
(owk 
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eidh 


where, [X] in Eqs. 4.1-4.3 is determined from the normalized 
Site conservation equation: 
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Since no simplifying assumptions have been made in deriving 
Eqs. 4.1-4.3, then the use of the site conservation equation 
to eliminate [X] is arbitrary, as any one of the four 
concentrations could be eliminated without incurring the 


effects described in the previous chapter. 


4.3 Stability Analysis 
The stability properties of Eqs. 4.1-4.3 are determined 


from the Jacobian matrix, which for this system 1S given by: 


-K,-K_,-4K,[AX][Bxx] -K, =2Ke=2KebAxi)2 
a= -K, -K,-1 -~2K, Us, 50 


Ree. WAX) PRxx | ke x ie ike | MRE se = 2h TARY ? 


The eigenvalues of Eqs. 4.1-4.3 are solutions of the charac- 


teristic equation associated with J, which is given by: 


—(eiwAGw@: —BwAsbidwhist=eo (236) 

where A,, A, and A; are the three invariants of J, 1.e.: 
Aye= Traced) er: 
Pope eM UN OTN ep at IME Or (ol oe ee Manors (3 Sa) (4.8) 


Roe= Determinant cd) 


andy thesterms, 748, 8 sheandejs6) ares theme! enentsmonsenecemarn 
diagonal of the Jacobian. From Routh's stability criterion 
[4.1], the necessary and sufficient conditions for a steady- 


state to be asymptotically stable are that: 
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Condiercnias | ae a0 (4.9) 
Condition #2 2A tet / Ae tO CARO} 
Cond ution ee je ae) C4544) 


If any of these conditions are violated, then the steady- 
State will be unstable, as at least one eigenvalue will 
have a positive real part. 

It 1S seen that Condition #1 is never violated since 
the Jacobian matrix only contains negative elements. Thus, 
the stability of Eqs. 4.1-4.3 will be solely dependent on 
Gong ie vomsy 2 sand (eee iwha le thembehnaviortof A ‘can be 
readily investigated analytically, no such investigation can 
be performed for A; due to the complexity of the terms in 
the Jacobian. Thus only a numerical stability analysis can 
be performed for this system as both Conditions #2 and #3 
are dependent on the value of Aj. 

A numerical stability investigation, which is not 
claimed to be exhaustive, has shown that when Condition £3 
1s violated at a steady-state, then two additional 
Steady-States exist for which Condition #3 1s not violated. 
The stability properties of these additional steady-states, 
found by examining Condition #2, determine the global stabi- 
lity icharactertstics of the system.” This results anti 
tively correct, as a steady-state is locally unbounded if 
A,>0, as this is the three dimensional analog to a saddle 
point. In a system obeying conservation balance 


constraints, an unbounded global solution is not possible, 


a 
ia 26 . GB >2.0A 
{ar “1 a ae : 
a , binds 
; giq@ela whoo Petal Sly G24 er-ov Th oer saeae’ to 138 
roiw qua eviegie a 3s pees nik wt iis He ; ie 
a | 7 
“ma Lob? subskeog 6 svhit Yh | 
| ee 
b ee ? 2 os yloe is Cnesetees ont , 
| J 


| . | b~t & . sell So Qe eee aaa 
| vi , ; J A al) ee ‘4 eis ty seols (Se f 
ity jee S00 eek ¢o eee 
+ al 4 Tete mee) 
| iy ‘teed Holts ‘es Co yRie wits iy Seeae? “ae 


ia oy ~ebaeasl EAP o 


a 
iz 


jvel *Mtidae Last oemee «& | ole 
ao 
- 


re 
= 


i + phoamug ee ae nauk ay 


+g e"¥hoole> &,. 8G hate ne? ale 
ayy: 


“i 


‘ 


: seeeeeoly Sit i SAT iaty hae The > ieive ewbusta 


Be unter ween ino: ete apie 476c7>y5 sy ge 


tivity Md yet 9 ih ty weletbes? £2 De: ea Oe 
" 7 


aes uae pet Reps it in > corde oF bei | 


on Kite re ay i _ +e ; : Ae 


. (Toe ¢ 


i 
{ 
0 
: 
% 
_ 
iy : 


hence if Condition #3 is violated at a steady-state, then 
additional steady-states must exist, as a limit cycle cannot 
Surround a unique steady-state which is a saddle point. 

In order to compare the general model composed of Eqs. 
4,1-4.3, with the simplified model presented by Eigenberger, 
ehempatvametercekKe Kw, Ke atoehke: will sbe held constant, 
where, 


Kee ore ok hee Ke ce eee 


Since these are the only parameters (in modified form) which 
are present in Eigenberger's simplified model, then any 
differences in behavior between the general model and the 
Simplified model will be solely due to the effects of the 
equilibrium assumption. These effects can be determined by 
examining the behavior of the general model in the K,-K_, 
parameter space, as these are the two parameters in the 
general model which result from the removal of the equili- 
brium assumption. 

A detailed numerical analysis of the stability proper- 
ties of the general model for the parameters K,=400, 
K_,=5.8309, K,,=20000, and K,=5 was performed. These values 
correspond to Eigenberger's [4.2] parameters k{,=0.02, 
[AX]g=0.7, ki=5 and Ta=20000, for which the simplified model 
behavior consists of a sustained oscillation about a single, 
unstable, steady-state, as shown in Figure 4.1. The 
behavior of the general model in the K.,.-K_, parameter Space 


is shown quantitatively in Figure 4.2 and described 
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Figure 4.1 
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Predictions of Eigenberger's Model 
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Figure 4.2 


Stability Regions for the General Model 
dSeagruncGtioneot Ko) andekK os. aae 
(Parameters as in Figure 4.1) 
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qualitatively in Table 4.1. The computer source code for 
the program which was used to produce Figure 4.2 is included 
in Appendix A, 


Tablevwe. 32 atTypes or Solutions: to Equations 4.1 = 4.3 
(Regions refer to Figure 4.1) 


Region Solution Type Stability Conditions 
A, TG NGA S A; 
1 one stable steady-state a = es 
ie one unstable steady-state = 3 = 


Giimitecycle)) 


three unstable steady-states 
(limit cycle surrounds 
all three steady-states) 


one stable, and two 
unstable steady-states 


In Table 4.1 it 1S seen that the general model is capable of 
exhibiting complex behavior not observed in the simplified 
model (for the previously indicated parameter values). In 
Pacvercllbar, cepending on the values of K, land Ki2, the 
general model can have either one or three steady-states 
with various different types of dynamical behavior. It is 
oe ae to note that this diversity of behavior as ail 
occurring for the situation in which the simplipiedsmodei s 
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4.4 Numerical Results and Discussion 

Numerical integrations of Eqs. 4.1-4.3 with various 
values of K, and K_, were done. Figure 4.3 shows the 
predicted behavior for values of K, and K_, located in 
region 2 on Figure 4.2, and Figure 4.4 corresponds to K, and 
K_, values located in region 3. 

A comparison of Figures 4.3 and 4.4 with Figure 4.1 
shows that when solutions of the general model display 
oscillatory behavior, these oscillations can vary signifi- 
cantly in period and amplitude from those of the simplified 
model. These differences are important, as one objective of 
modelling 1s to obtain a model that quantitatively, or at 
least semi-quantitatively predicts experimentally observed 
behavior. Since the experimental data that are most easily 
Sbeamneceaclemriewintinences.of Operating conditions on) ehe 
period and amplitude of the oscillations, then the first 
test of the adequacy of a model is whether or not it 
correctly predicts these effects. Therefore, the frequency 
and amplitude predictions of a model are very important in 
testing the validity of the model. 

Perhaps aS important, Figure 4.2 shows that the general 
model displays stable behavior and multiple steady-state 
behavior for large ranges of parameter values over which the 
Simplified model predicts oscillatory behavior. Hence, the 
general model predicts behavior which is qualitatively and 
quantitatively different from that displayed by the 


Simplified model for the same parameter values. This 


te 


er we 


Mh as ead ees? abe Gow 3 


Bip siniw ae 


pi) Oa He | 
A ugt’ ne & oo ion’ 
~~) Q—gaeet sraler eis 
a not gs = 
ulev nage todd aga 
vases 2 ul i "San 


ner ¢ age) Se '. Jae 


;, 4 ~)4,h 60@** | «590Rm 
obs (ooee 
iieup~ inet. Seems 


Sereer fi ’ iff 4 wa | _ al 


ae 
Si we 16 Geol ste _ 
ne ae 
fi) to souorlqne Pre boise " 7 
ts » csaepels ectd Do) SRS 


sane: si 3ibemty. 2i23 2sA169) oe 


» in oro) esi bet agate onl Cae 


ener Jacds. Rae pt tke ee oie 


it 
‘ort, (ie heen oe squeeey 
ye acne ang ie cere 


ie LP 


ro 


Li 


[AX] 


[BX] 


[x] 


Fractional Surface Coverages 


[BXX] 


Figure 4.3 


66 


Dimensionless Time, T 


Predictions of the General Model 
(K,=6, K.,=10‘*, other parameters as in Fig. 4.1) 
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Dimensionless Time, T 


Predictions of the General Model 
(K,=1, K.,=10‘, other parameters as in Fig. 4.1) 
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indicates that assumptions regarding equilibrium behavior of 
certain reaction steps can greatly alter the predicted 


behavior of a system. 


4.5 Conclusions 

It has been demonstrated, contrary to the claim of 
Pereira and Varama [4.3], that the use of the equilibrium 
assumption in the construction of a reaction model can 
qualitatively and quantitatively alter the predicted 
behavior of the system. This particular simplifying assump- 
tion can affect both the number, and the stability, of the 
Steady-states for a reaction model. 

Great care must be taken when assumptions are intro- 
duced into a model to ensure that these assumptions do not 


obscure the behavior of the original model. 


4.6 Nomenclature 


Roman 
A; .coefficients of the characteristic equation 
[A] gas phase concentration of A 


[AX] normalized surface phase concentration of A 

[B] gas phase concentration of B 

[BX] normalized surface concentration of reactive B 
[BXX] normalized surface concentration of unreactive B 


is element on the main diagonal of the Jacobian matrix 
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Jacobian matrix 

dimensionless rate constants for the forward reactions 
dimensionless rate constants for the reverse reactions 
combined equilibrium and reaction rate constant 
dimensionless time 


normalized active site surface concentration 


eigenvalue of the characteristic equation 
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5. CHAOS AND RELATED PHENOMENA 


Sein eroduc ts on 

In the preceeding chapter, the behavior of a general 
model consisting of three ordinary differential equations 
was compared to that of a simplified model consisting of two 
ordinary differential equations. It was seen that the 
steady-state and dynamic behavior of the general model was 
considerably more complex than that of the simplified model. 
However, due to the ease of analysis, most of the recent 
work in examining the stability of chemical reactors has 
been directed at reactor models which have consisted of two 
ordinary differential equations. The variables considered 
in these models have usually been either two reactant 
eoncentrations 1521) 5.2], or one reactant concentration and 
thewreacton temperattrem( Sas =£:56°5) . 

One of the implications of examining systems of two 
ordinary differential equations, such as Bigenberger's 
Simplified model, is that the most complex dynamic behavior 
that can be obtained is a sustained oscillation of the two 
variables, which can be represented as a limit cycle in the 
phase plane. However, recent experimental studies [5.6 - 
5.8] have shown that chemical reactors Can exhibit’ chaotic 


behavior (cycles with high periods, or even aperiodic 


») 


behavior [5.9]) under certain operating conditions. I 


_ 


order to model such behavior, it 1S necessary to use a 
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minimum of three ordinary differential equations, as the 
Poincare-Bendixson theorem prohibits the existence of 
chaotic behavior in systems of two ordinary differential 
equations. This theorem states that a bounded solution of a 
two-dimensional dynamical system will converge to either a 
Stable steady-state, or to a stable limit cycle [5.10]. 

The modelling of chaotic behavior in chemical reaction 
Systems VSsistill in Gis infancy f as only a few studies (5.8, 
ol pee oeemeinave sShown Bthat realrstic weaction models can 
Gasplaye chacticwso lutions gel n ‘the ifol lowing “it twin ibe 
shown that chaotic behavior can be exhibited by a very 
Simple chemical reactor model which consists of three 


ordinary differential equations. 


5.2 Model Development 
The parallel conversion of two reactants A and B into 
respective products C and D can be written as: 
A —> C eee By, 
B —— OD Cae) 
If these reactions are carried out in a nonisothermal, 
continuous flow, Stirred tank reacton, thenwtnis chemical 
system can be modelled by a system of three ordinary 
differential equations which result from the application of 
the principles of mass and energy conservation. In deriving 
the model, the following assumptions will be used: 
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2. The reactor is well mixed so that the temperature and 
the species concentrations are uniform throughout the 
reactor. 

3. Both reactions are exothermic, and follow first order 
rate expressions whose temperature dependencies obey 
the Arrhenius equation. 

4, There are no volume changes on reaction. 

5. The physical properties of the reaction mixture are 
temperature invariant. 

Gam thetreactiron fluid phase contains the only significant 
thermal capacitance. 

7. The temperature of the reactor wall is held constant at 
the feed stream inlet temperature. 


GeetTie reactor wall heat transfer coefficient is constant. 


Using the above assumptions, reactor mass balances for 
Species A and B, and an overall reactor energy balance, take 


the form: 


-—,/RT' 
Vio kAd es OC [AW (Ad)) —2 VLA kee 58S) 
sian 
-E,/RT' 
v a[B] = O([B].-[B]) - V[IB]k.e (5.4) 
ots 
-E,/RT' -E,/RT' 
VeCp dT' = Hr,V[A]k,e + Hiv. Lb ese 
ot 
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The three balance equations, Eqs. 5.3-5.5, can be written in 


the following dimensionless form: 


vel, ale T) 
dva=wriee Yoruypaye coo 
dt 
ree AC ae 
Clg Mess cael Daye on) 
at 
Veet y,T/(1+T 
aT (B,YDa,e Pap y2Da,e mu eee ch) 
dt 
Liv 
where, Da, = Vk.e (oa) 
Q 
Sate 
Da, = Vk.e (5.9b) 
Q 
Y¥4 =_Ey Cee) 
RT! 
V2 = BE, (5.9d) 
RT} 
bo ="O\Hry TA]. (5.9e) 
Tt (SU+OpCp) 
B, = Q(Hr,){B]\, (5.9£) 
wis (Si+0eCp) 
ie - slext eeu (S.0ce 
QpCp 
ta. Om b5 9h) 
V 
Fe eee CBS 1) 
ee 
y = [Ad (5.295) 
eat. 
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An identical set of equations, with different physical 
interpretations of the parameters, is obtained for two 
parallel reactions occurring in a nonisothermal, heterogene- 
ously catalyzed system if the surface concentrations are 


proportional to the gas phase concentrations. 


5.3 Stability Analysis 
The stability properties of Eqs. 5.6-5.8 are determined 
from the Jacobian matrix, which for this system of equations 


1S: 


eT Aas) Stil ioe tan e | 


-1=Da,e 0 Ly Dare 
(T+1)? 


eet at ily) velar a Slay) ee 
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As previously shown, the eigenvalues of Eqs. 5.6-5.8 are 
solutions of the characteristic equation associated with J, 
which is given by: 

SAR oe Ba Oey mine Ae 10) CBee a) 
where A,, A, and A, are the three invariants of J. The 


necessary and sufficient conditions for a steady-state to be 
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If any of these three conditions are violated, then the 
steady-state will be unstable, as at least one eigenvalue 
will have a positive real part. 

A general analytical investigation of when the above 
Stability conditions are violated is not possible due to the 
complexity of the terms in the Jacobian matrix, and due to 
the possibility of multiple steady-state solutions of Eqs. 
5.6-5.8. A combined analytical and numerical uniqueness 
analysis [5.13] has shown that, for certain parameter 
values, aS many as five steady-state solutions can exist, 
each with potentially different stability properties. 
Certainly a dynamical model which possesses a multitude of 
coexisting steady-states with a diversity of stability 
properties could be expected to exhibit complex global 
behavior. However, we shall show that when the parameter 
values are restricted so that the steady-state is 
necessarily unique, both regular and chaotic cycling may 


still ensue. 


Sree Hopf Bifurcation from a Unique Steady State 

Using the techniques outlined by Luss [5.13], it is 
possible to show that Eqs. 5.6-5.8 will have a unique 
steady-state over wide ranges of parameter values. In 
Darticular, if Dag=2. 7/5) 6, =0.0%,) y.-92-2o ano o=2o0,. Uren 


a unique steady-state exists for all values of Da, and £,. 
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For the above parameter values, the stability and 
character (node, focus, etc.) of the unique steady-state in 
the Da,-B, parameter space can be determined by calculating 
the eigenvalues of the characteristic equation. Shown in 
Figure 5,1 are the regions of equivalent character of the 
unique Staedy-state, where it is seen that the steady-state 
can be one of four types: (1) an unstable focus (one nega- 
tive real eigenvalue and two complex conjugate eigenvalues 
with positive real parts); (2) an unstable node (one nega- 
tive real eigenvalue and two positive real eigenvalues); 
(3) a stable focus (one negative real eigenvalue and two 
complex conjugate eigenvalues with negative real parts); or 
(4) a stable node (three negative real eigenvalues). 

The manner in which a stable steady-state becomes 
unstable as the parameters are moved from one region to 
Shonncte swShOWn ine rigure 6.2. "in thrs figure, the para- 
meter B, 1S moved from region 1 to region 3, and the effect 
of this on the time behavior of the dimensionless tempera- 
ture is illustrated. An examination of the eigenvalues 
listed in Figure 5.2 reveals that a pair of complex conjug- 
ate, purely imaginary eigenvalues occurs at the boundary 
separating region 1 from region 3. At this boundary, it has 
been numerically determined that the real parts of the 
complex conjugate eigenvalues have non-zero derivatives with 
mespect «tO B... Thus) .a Hopt .bi furcations | ow eaoccurs of 
the boundary separating region 1 from region 3. An 


examination (not shown here) of the eigenvalues as £, is 
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Figure 5.1 Stability Regions as a Function of Da, and B, 
(Da, = 2.75, B, = 0.040, a = 250, y, = y2 = 25) 
Insert shows location of parameters for Figures 5.2 to 5.7 
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moved from region 2 to region 3 reveals that a Hopf bifurca- 
tion also occurs on this boundary. 

ineeallecases 1 he numerical sntegration of Eqs. 5,675.8 
was performed by using Gear's [5.15] method for the integra- 
tion of stiff systems of ordinary differential equations. 
For most cases, the results were also checked by uSing a 
fourth order Runge-Kutta-Fehlberg [5.16] method, with fifth 
order error estimation, to integrate Eqs. 5.6-5.8. Hence, 
two very dissimilar integration methods were used to 
integrate the system of equations, and similar results were 
always obtained. The computer source code for these 


programs is included in Appendix B. 


5.5 Complex Oscillations and Chaotic Behavior 

From the /work of Lorenz [5.17], it is now well known 
that systems of three or more differential equations can 
exhibit solutions of arbitrarily high period, as well as 
aperiodic solutions which are not even asymptotically 
periodic. This type of behavior is frequently referred to 
as chaotic behavior. For Eqs. 5.6-5.8, parameter values 
have been found for which the system can display apparently 
chaotic behavior. The transition from complex periodic 
behavior to chaotic behavior, and then back to complex 
periodic behavior 15 ishown }in Figuress oe comonwe. 

Bach Ob Figures os to 5.6Mare omposedmor the 


following three sections: (a) a plot of the time behavior of 
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the dimensionless temperature, (b) an autocorrelation 
analysis of the temperature maxima, and (c) a next-amplitude 
plot for the temperature maxima. In the following, each of 
Pune. Guoupseofs (a). ib) eandy(c)ssectionsefon Bigures:5<3:to 
S.6pwill be discussed/ineturn, 

DnwBagquress5.3(a)eten5s6(a)) thesvariation -ofsdimen- 
Slionless temperature with time is shown. The behavior for 
the first 26 units of dimensionless time has been discarded 
cCoreliminate thestransient effects of the initial condi- 
BrOns Ml neeccnecace the inibraly conditions are T3=+0, Y,=0, 
and Z,=0. We see that, as $B, 1S increased, the simple 
Oscrllaticne displayed) in Figure 52 brturcates into the four 
Naw oeperecvc Lemoscimetron of Fiore 5,344). “AS pp, “is 
further increased the temperature is shown to behave in a 
chaotic fashion as in, Figures 5.4(a) and 5.5(a). An addi- 
tional increase in 8, results in the temperature once again 
following a stable periodic trajectory, this time with six 
maxima per cycle, as shown in Figure 5.6(a). Besides stable 
periodic oscillations of orders four and six, it is possible 
to have many other numbers of maxima per cycle (2, 3, 6, 8, 
9, 16, etc.) depending on the value of B.. In many systems 
which display chaotic behavior, it has been found [5.18] 
that period doubling occurs as the system behavior changes 
Epomea Stable periodic to a chaotic solution. 9 Sinwgeneral, 
for the system of Eqs. 5.6-5.8, period doubling does not 
occur, although it has been observed for some sets of 


parameter values. 
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Mme Firqures 5.3 ¢b)ete 5.6(0b) the*results*of an auto- 
correlation analysis on the temperature maxima are shown. 
The autocorrelation analysis was performed by first reducing 
the complete time series of temperatures obtained from the 
integration program, into a vector of temperature maxima, 

Ti Gmeemwienemi Varteoerrome tetoOlN, Vand @Nls4the®number of 
temperature maxima in the time aaron The autocorrelation 
coefficient for L lag units (correlation between Tm(i) and 


[merce ewas chemecalcu lated "from: 


N-L 
AG (ih) e-em eee em (Tm) Tm )4tm Cit) Tm). (oi) 
NST. Var (Tm) 
1=1 
ne N 
where, nmes ee etm () (592) 
n / 
i=1 
N _—— 
Vou (Pi) vie erm (a) Tm)? CSe4) 
NH / 


It is important to realize that the temperature maxima are 
considered to be equally separated by one lag unit when they 
are stored in the vector which only contains the temperature 
maxima, even though the maxima were perhaps separated by 
varying amounts of dimensionless time in the original 


unreduced time series. 
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In Figure 5.3(b) the autocorrelation analysis clearly 
indicates that an oscillation with four maxima per cycle 
occurs because the autocorrelation coefficient is equal to 
iciueasVa LlCaGreboursbaq unitseand at allimultiples of 
DOUREUOGMUNI tS wi thewlatrgemvarlve of the autocorrelation 
coefficient at a value of two lag units indicates that the 
four maxima per cycle oscillation is very close to being a 
two maxima per cycle oscillation. \.This type of behavior is 
also shown in Figure 5.6(b) where the six maxima per cycle 
oscillation is clearly shown by the autocorrelation coeffi- 
Cle temavi Onda VaMuerOred eat) SLXalaq unitSssandjatuall 
Multrplessot.si x laqwunits. 

For the chaotic behavior displayed in Figures 5.4(a) 
Sndeon5(a)yjescels seen that the corresponding-auto- 
correlation coefficients given in Figures 5.4(b) and 5.5(b) 
are never equal to one, and in fact they approach zero at 
large lag values. While only twenty-four lag units are 
SHOWNeIn Pigures.5.4(b) sand 5.5(b) , ethesautocorrelation 
coefficients have been calculated for several hundred lag 
units, and no long periods of oscillation were detected. 
This particular form of autocorrelation analysis, as 
detailed above, has been found to be a very useful, and 
computationally efficient, method of detecting chaotic 
Hbehavuors in the aso lucvon otha Ss 60 .6> omc. 

Pnakiouresms. 3 (Cc) vtGno.6 (Cl sale wSnOWnmunoume xi 
amplitude plots, in which successive dimensionless 


temperature maxima Tm(it+1) are plotted against Tm(i). If a 
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Stable periodic solution exists, this will be shown on a 
next-amplitude plot by the existence of a finite number of 
points which correspond to the number of maxima per cycle in 
Biesperog1c Oscillation. for example, Figure 5.3(c) has 
COMME LOUr pOIntS ONelt ana this corresponds to an oscilla- 
tion with four maxima per cycle. Similarily, the six maxima 
per cycle oscillation in Figure 5.6(a) generates a next- 
amplitude plot containing exactly six points, as shown in 
Picures 536.0 ):, 

The next-amplitude plots become more complex for the 
case of chaotic oscillations. Lorenz, ina series of papers 
(cf. references in [10]) modelled turbulent fluid flow 
through a simplified set of equations. From these he numer- 
ically solved for X(n), the maximum kinetic energy of 
Successive waves. He plotted X(n+1) versus X(n), and found 
that a continuous cusp-shaped curve y=F(x) (with a single 
maximum value at the cusp) approximated this data. However, 
asecans be.seen sn higures 5.4(c) tand S25\c)}, the chaotic 
behavior of Eqs. 5.6-5.8 does not generate a Single curve, 
rather, at least two curves are visually apparent. This is 
not to imply that the relationship between Tm(i) and Tm(i+1) 
is non-unique, but rather, as Tm(i) varies there are an 
infinite number of switchings between the curves. The 
curves must be discontinuous because the next-amplitude map 
'S Single walued;. ive. a verticai line Whvehecrosses through 
a point on one curve necessarily crosses through no point on 


the other curve. Figures 5.3(c) to 5.6(c) do not represent 
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the full extent of possible forms of the next-amplitude 
plots, aS even more complex forms can exist depending on the 


choice of parameters. 


5.6 Chaos and Fractal Curves 

ShOwminehigureso, (eC) iS a nextramplitude plot in 
which three curves are distinctly observable. The chaotic 
oscillation of the dimensionless temperature which generated 
this next-amplitude plot is shown in Figure 5.7(a), and the 
autocorrelation analysis is given in Figure 5.7(b). 

Sections of the curves shown on Figure 5.7(c) were 
expanded to more closely examine their structure. Shown in 
Figure 5.8 are successive magnifications of the regions 
enclosed in the dashed rectangular boxes. Figure 5.8(b) 
shows quite distinctly that the "single" curve in the dashed 
box of Figure 5.8(a) is actually composed of at least three 
GCurvese 9 in Figures 5,8(c) and 5.8(d) at is seen thet upon 
additional magnification, these curves expand into even more 
Gurves. lt 1s reasonable to expect that “cunves son the 
next-amplitude plot shown in Figure 5.7 are actually 
composed of an infinite number of curves, and in fact these 
are fractal curves in the sense of being self-similar 
[5.19]. Specifically, the transverse structure across the 
roughly parallel curve segments is essentially that of a 
Cantor set, while on the curves the points appear to distri- 


bute themselves densely. 
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Figure 5.8 Fractal Curves in the Next-Amplitude Plot 
(a) Next-amplitude plot from Figure 5.7(c) 
(b) 278:1 expansion of region in dashed box in (a) 
(c) 9.4:1 expansion of region in dashed box in (b) 
(d) 8:1 expansion of region in dashed box in (c) 
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Such a topological structure has been encountered in a 
number of examples of second-order recurrences which 
approximate certain three-dimensional dynamical systems that 
have been used to describe turbulent flow. These recur- 
rences take the form of the Poincaré return map to the flow, 
i.e. the map which describes the flow on some transverse 
planar section. A well-known example is the Hénon return 
map to the Lorenz equations [5.20]. 

It should be mentioned that the production of Figure 
SeCagwasunotwantnhivialitask,rasi the approximatelye 1507 points 
shown in Figure 5.8(d) were obtained by integrating Eqs. 
5.6-5.8 so that approximately 15000 temperature maxima were 
recorded. This integration required slightly more than ten 
days of CPU time on a dedicated Hewlett-Packard 1000 21MXE 


minicomputer. 


Safer rajectories sin iT—-Y—Z2 Space 

While the manner in which the variables T, Y and Z 
change with time is of interest, the forms which they create 
in three dimensional T-Y-Z space are of even more interest, 
as these will delimit the shape of the attractor which the 
trajectories are following. For the chaotic oscillation 
ShoWwiean Eroure 5.7) the T-Y 1-24 andee vevilewscucusuue 
system trajectory are shown on Figures 5.9(a)-(c) respec- 
tively. From these two-dimensional views, the complete 


shape of the three-dimensional trajectory can be 
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constructed. 

Several positions have been labelled in Figure 5.9(b) 
in an attempt to describe the three-dimensional shape of the 
berecloOlmwe Ul saute ecrOony starts nNéar position | °on ‘the 
front (out of the paper) of the lower band, it will travel 
Prone Gone con selt. andmrascetnroughnepostTt rom. such’ that<it 
is now located on the back (behind the paper) of the upper 
band at position 3. From there, the trajectory passes under 
ENewr1ont mostepamibaotethewatbuactonssuch-that-ltiwarrives at 
position 4 which is on the back of the lower band. The 
trajectory passes along the inside of the "circular-well" 
(clearly shown in Figure 5.9(a)) in a clockwise fashion 
until it reaches position 57 which "1s "onsthe front side of 
the upper band. From position 5, the trajectory travels 
OvemeUnesCOUECLA Une Tr EGhtLomoste partworethe attractor, ultim- 
Seehveanmrving abaposutlon, 6, swiich Vike tposition i> 41s 
located on the front of the lower band. The trajectory will 
now repeat itself in a similar way. 

This motion is identical to that which would result 
from following a Mobius-band, hence this attractor should be 
known as a "Mobius-band attractor". In some respects, this 
attractor is similar to an idealized attractor described by 
Rossler [5.21]. Besides this particular shape, many other 
forms of the attractor have been observed for other 


parameter values. 
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Figure 5.9 Two-Dimensional Views of the Three-Dimensional 
Attractor (all parameters as in Figure 5.7) 
(a) T-Y view (clockwise motion when viewed from above) 
(b) T-Z view (Mobius-band structure is displayed) 
(c) Z-Y view (motion is clockwise) 
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Bee Conclusions 


It has been shown that a very simple chemical reactor 
model consisting of three ordinary differential equations 
Galec splaymcomplexs periodic and chaotic behavior. “This 
behavior has been analyzed by examining autocorrelation 
coefficients, next-amplitude plots, and state-space trajec- 
tories. 

Through the use of the autocorrelation analysis it has 
been shown that the chaotic behavior does not contain any 
high order oscillations. When only the temperature maxima 
are used, the autocorrelation analysis is a very useful, and 
computationally efficient, method of examining complex 
perlodrec®andcchaotic* behavior. 

The next-amplitude plots have been found to be composed 
of very complex structures, which upon close examination 
appear to be fractal curves. 

A detailed inspection of the state-space trajectories 
auringtachactic*oscillation has shown that the trajectories 
follow a Mobius-band attractor. Other more complex forms of 


the attractor have also been found. 


5.9 Nomenclature 
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Greek 

a dimensionless heat transfer coefficient 
Be dimensionless heat of reaction 

Yi dimensionless activation energies 

p density 

Ww eigenvalue of the Jacobian matrix 
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6. DESCRIPTION OF THE EXPERIMENTAL APPARATUS 


In the preceeding chapters several topics have been 
investigated which are of general interest in the modelling 
of oscillatory behavior in chemical and catalytic reactors. 
However, no attempt has been made to exhaustively 
investigate the large number of recent models which have 
been proposed as explanations for oscillatory behavior. It 
appears that the devising of mechanisms for explaining 
oscillatory behavior is only limited by the imagination of 
those who model. In order to determine which mechanisms are 
physically realistic, it is necessary to perform model vali- 
dation through comparison with actual experimental observa- 
tions. To this end, a reactor for studying the dynamics of 
chemical reactions over supported metal catalysts was 
constructed. This reactor and associated equipment are 


descrabedsin this chapter. 


6.1 General Description of the Equipment 

As initially conceived, the experimental apparatus was 
to consist of a recycle reactor system which could be used 
to study the oxidation of carbon monoxide on a supported 
platinum catalyst, as this was one of the reactions for 
which there were reports of experimental observations of 


oscillatory behavior. Since previous experimental work upon 
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which to base a design was not available, it was necessary 
to adhere to several very general considerations in the 
design and construction of the apparatus. 

It was desired that the apparatus be capable of being 
used for general investigations of gas phase reactions on 
Supported metal catalysts, i.e. 1t was not to be constructed 
specifically for the oxidation of carbon monoxide on one 
particular catalyst. In particular, temperature measurement 
and control over a wide range of operating conditions was to 
be possible, and similarily, the gas metering and delivery 
system was to be capable of accurately providing gas 
mixtures over a wide range of compositions and flow rates. 
Also, due to the dynamic nature of the phenomena which were 
to be studied, continuous product analysis was essential. A 
general description of the apparatus which was constructed 
PSeGuvenwitiFrgqure sao, ancsin the oltowing, seach of the 


major sections of the apparatus will be described in detail. 


6.2 Feed Metering and Mixing 

Shown in Figure 6.2 is a detailed schematic of the feed 
ioe he and mixing section of the apparatus. As can be 
seen, provision was made for using various gases (0.,, CO, 
NN) He, é6tc..) for the purposes of reaction or dilution. 
With one exception, all these gases were Supplied from high 
pressure cylinders obtained from Union Carbide Limited 
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Figure 6.1 General Description of the Experimental Apparatus 
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Limited. The oxygen and nitrogen were very dry grade (99.6% 
and 99.99% min. purity, respectively), the helium was 
prepurified grade (99.995% min. purity), and the carbon 
monoxide was: C.P .qrade (99.5% min. purity)... * The: outlet 
pressures from the gas cylinders were controlled through the 
use of several Matheson model 8H two stage regulators. The 
sole exception to the above was a aves 38 mlitnas huch 
purity oxygen (99..99% min. pouraty)\ obtained from the 
Matheson Company, the pressure of which was controlled via a 
Matheson model 19 high purity, single stage, regulator. 

The flow rates of the various gases were measured by 
two Matheson rotameters (#600 and #610 tubes) and by two 
Matheson linear mass flowmeters, model 8116-0153 (Ser. 
fie FO) e.0-5000.SCCM, andamodéels ALLA 005.(Ser.4 241002), 0-106 
SCCM. The voltages generated by the two mass flowmeters 
were continuously recorded by a two pen Hewlett-Packard 
model 7100B-14-21-25 (Ser. #844A04091) chart recorder. The 
flow rates of the gases were adjusted by using metering 
valves located at the exhaust end of the rotameters, and by 
two Whitey model 21 fine metering valves located on the 
downstream side of the linear mass flowmeters. The pressure 
in each rotameter was constant at the outlet pressure of the 
cylinder regulator which supplied the gas since all the 
valves were located downstream of the rotameters. 

The flowmeters and rotameters were calibrated for a 
number of gases at two different pressures, and these 


calibrations are included in Appendix C. The linear mass 
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flowmeters were found to be insensitive to pressure over the 
range of the calibrations and the high range mass flowmeter 

(0-5000 SCCM) was only calibrated for a small range of flow 

rates since a 50 cm? bubble flowmeter was used for the cali- 
brations. 

From Figure 6.2 it is seen that if the gases were not 
sent to the bubble flowmeter, then they were instead passed 
through a particulate filter (particles seven microns in 
diameter or greater were removed) and a static gas mixer. 
The gas mixer consisted of stainless steel baffles placed 
inside a 20 cm piece of 1.9 cm diameter stainless steel 
tubing. From the mixer, a three way ball valve was used to 
direct the gas either to the recycle compartment, and on to 
the reactor, or to the gas chromatograph and the infrared 
spectrophotometer for analysis or calibration. In the 
construction of the feed metering and mixing apparatus, 
brass tubing and valves were used upstream of the rotameters 
and flowmeters, and stainless steel valves and tubing were 


used on the downstream side. 


Gao nested Recycle Compartment 

Shown in Figure 6.3 1s a detailed schematic of the 
heated recycle compartment. The primary purpose of the 
recycle compartment was to provide feed to the reactor at 
approximately the reaction temperature. The temperature in 


this compartment was measured with an iron - constantan 
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thermocouple (type J), controlled with a Honeywell model 
R7161J (Ser. #1015-1) temperature controller, and recorded 
on a Honeywell model 687378-1 (Ser. #68-06-07) Electronik 16 
MuUPDI point we sAuthinty-three ohmiresistance coil wath a 
radiation shield was used as the load for the temperature 
controller, and the heated air was circulated by a Universal 
Electric oupmodel 80-0419 (Ser. 46A54775R) motor with a 
Squirrel cage fan. Heat losses from the recycle compartment 
were minimized by placing fiberglass insulation on the 
inside walls of the compartment. 

Located in the recycle compartment were the necessary 
pieces of equipment for recycling the reactor product stream 
and controlling the flow of the effluent stream. Stainless 
steel was used throughout as the material of construction. 

The recycle was accomplished with a Metal Bellows 
Corporation model MB-118HT (Ser. #1622) metal bellows pump. 
The metal bellows pump was driven with a variable speed one 
horsepower Reliance VSD motor, model B56G3104 (Ser. 
#27F85W05-AZ). This motor was manually controllable over 
the range 30 to 2000 RPM by means of an associated 
electronic speed controller. A gas filter capable of 
removing particulate matter with diameters of seven microns 
or greater was placed on the suction side of the recycle 
pump in order to protect the metal bellows. As indicated in 
Figure 6.3, the flow rate of the recycle Stream was moni- 
tored with a rotameter (Matheson #605 tube), and the 


pressure in the recycle line was measured with an Ashcroft 
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(30 inches Hg vacuum - 30 psig) pressure gauge. 

The flow of the effluent stream was controlled with a 
Moore constant differential flow controller, model 63SUI 
(Ser ePG- 5402), # Thiss flow conéroliléer was used to 
manipulate the effluent flow rate such that the pressure in 
the recycle line (and in the reactor) was kept equal to an 
externally applied pressure. This external pressure (set 
point for the flow controller) was provided by a Moore 
Nulmatic model 40-100 (Ser. #PC-7170) pressure regulator 
which was connected to a low pressure air line. 

While a gas analysis was normally performed on the 
reactor effluent stream, it was also possible, as shown in 
Figure 6.3, to sample the stream (recycle plus new feed) 
entering the reactor. However, sampling this stream was 
Giprr cult las) Wawas not sunder the control of Lhe flow 
controller. This line was primarily used to sample the 
reactor effluent when the reactor was run aS an integral bed 
reactor. For integral bed operation, the recycle pump was 
turned off, and the flow direction was opposite to that 


indicated in Figure 6.3. 


6.4 Catalytic Reactor in Heated Fluidized Sandbath 

A schematic of the fluidized sandbath, and catalytic 
reactor in which the catalyst particles were placed, is 
shown in Figure 6.4. The temperature of the inlet feed to 


the reactor was controlled by a Thermo Electric 400 model 
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Figure 6.4 Reactor in Sandbath 


eRe AEE 


iennaucse 4ee A eke io = 


=~ > Nietelieiameanil ee -_ - . B 
t _ pom ~ ntl ee il g 4 
int to bs 

7 ' 70 + ad 
‘ = ° 


=) 
a 
& 
= 
= 
= 


esiav 


f 
= 
Pay e> 
4 
aoe 
aes" 
Mpa bint 9) 
= y 
ial iat - 
a 
a 
aes 


ee 


>| — ey a at 
5 
ae 3 z 
- 
a = sand 
7 ae - - a 
~ ‘ 
- ¥ r 
ot “ees 
i; i 


108 


3242200 (Ser. #56486-118) temperature controller which was 
used to heat the sandbath walls. These walls were wrapped 
with resistance wire which presented a fourteen ohm load to 
the temperature controller and the resistance wire was 
covered with insulation to minimize heat losses. The heat 
released at the sandbath walls was transferred to the 
reactor by the fluidized sandbath. The sand (alundum) was 
fluidized by passing air through a sintered metal plate at 
the bottom of the bed of sand and a fine mesh screen was 
clamped to the top of the sandbath to prevent excessive 
losses of sand. 

The details of the reactor construction are shown in 
Figure 6.5 where the placements of the catalyst bed and the 
various thermocouples are shown. Stainless steel was used 
throughout as jtthe material of construction. With one excep- 
tion, ell of the thermocouples, including those welded to 
the reactor wall, were type J (iron - constantan). The sole 
exception was a type T (copper - constantan) thermocouple 
located in the tubing at the end of the reactor preheat 
section. | Thais? location’ was he control point for che 
temperature controller which required a type T thermocouple, 
hence the different type. The voltages produced by the 
thermocouples were recorded by the previously mentioned 
Honeywell multipoint, and in addition, the voltages of the 
control thermocouple and the thermocouple located 
immediately below the catalyst bed were continuously 


recorded using a two pen Varian model A-25 (Ser. #49044924) 
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Figure 6.5 Reactor Detail 
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chart recorder. 

The total volume of the empty reactor, the recycle pump 
(operating at 1800 RPM), and the associated tubing was 
determined to be in the range 245 to 260 cm*. This range of 
volumes was determined by observing the effluent carbon 
monoxide concentrations from a step input of carbon monoxide 
in the feed to the empty (no catalyst) reactor. The volume 
of the reactor was then calculated by assuming that the 
reactor and tubing could be approximately described as a 
first order system with a time delay. Due to the size of 
the reactor, a maximum of approximately thirty grams of 
Supported (Al,0,) catalyst could be used. This catalyst was 
Supported by a fine screen at the bottom, and retained by a 
balleote glass) wool at the top. The flowsdirection in the 
reactor was downwards when operated as a recycle reactor, 


and upwards when operated as an integral bed reactor. 


6.5 Auxiliary Reactor 

The auxiliary reactor, shown schematically in Figure 
6.6, was designed and constructed in order to improve the 
accuracy of the measurement of the concentration of carbon 
monoxide in the reactor effluent. This was achieved by 
completely converting all of the carbon monoxide to carbon 
dioxide. The concentration of carbon monoxide was then 
determined by taking the difference between the carbon 


dioxide concentrations in the effluents from the auxiliary 
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Eéeactor ‘and’the main veactor. “This resulted in an improve- 
ment in the accuracy of the measurement of the carbon 
dioxide concentration due to the infrared spectrophotometer 
being an order of magnitude more sensitive to carbon dioxide 
than to carbon monoxide (as will be shown in the next 
section). 

The construction details of the auxiliary reactor are 
quite straightforward, as the reactor was made from a 25 cm 
Piece of 1.2 cm diameter stainless steel tubing. Approxim- 
ately ten grams of 4 wt% Pt/Al,0O, catalyst was placed in the 
reactor, and the remaining volume was filled with glass 
beads. Due to the elimination of most of the dead volume 
through the use of the glass beads, it was found that the 
auxiliary reactor was also of use in determining some of the 
features of the dynamic behavior of the main reactor. 
Further details with respect to this are provided in the 
following chapter. Two gas filters capable of removing 
particulates with diameters of seven microns or greater were 
placed at either end of the reactor to prevent the contami- 
Matvon., col Ghe connecting tubing withetine catalyst 
particles. The gas from the main reactor was routed to the 
auxiliary reactor by opening one block valve and by 
Switching one three way valve. All valves, filters and 
bUDINGg were constructed OL (stainless Ss recis 

The temperature in the auxiliary reactor was measured 
by a type J thermocouple placed in the catalyst bed. This 


temperature was controlled by a Honeywell model R7161d (Ser. 
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#1007-1) temperature controller, and it was recorded by 
using the Honeywell multipoint. The load for the tempera- 
PULGEcOCneroblern consisted ef viortyyzone ohms of “high resis- 
tance wire which was tightly wound around asbestos 
insulating tape on the reactor. For safety purposes, a 
thick layer of insulation was placed over the resistance 


Waite. 


6.6 Gas Composition Analysis 

The analysis of the compositions of the various streams 
was accomplished by using both a gas chromatograph (GC) and 
an infrared spectrophotometer (IR). The arrangement of 
valves and tubing (stainless steel throughout) connecting 
the GC and the IR to the rest of the apparatus was such that 
any gas stream could be completely or partially vented, or 
could be sent to the GC and/or the IR. In the experimental 
work that was performed, the IR was used as the main tool 
for collecting dynamic information due to the continuous 
nature of its output, and the GC was used primarily for 
determining the overall compositions on a discrete basis. 
The use of the GC was necessary Since Several of the 
components in the reactant and product streams (0., N. and 
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6.6.1 Gas Chromatograph 

AS part of the experimental equipment, a Hewlett- 
Packard model 5734A (Ser. #1825A01203) gas chromatograph 
was used. This GC contained both a flame ionization dete- 
Secon raidwas thermal sconductivipy =a(TO)wdetector. ln thisuwork 
only the TC detector was used. 

The gas sampling was performed through the use of a gas 
sampling valve which contained two sample loops (0.5 cm® per 
loop). The sample loops were arranged such that one loop 
was continuously purged with sample when the other loop was 
used to inject a sample. Thus, there was a continuous flow 
of sample through the sampling valve. 

The output from the GC was recorded and integrated by a 
Hewlett-Packard model 3380A (Ser. #1324A01206) electronic 
integrator. This integrator was very user oriented, as with 
the push of one button it would automatically (under the 
control of the various switch settings) record and integrate 
the GC signal, produce a printed analysis of the signal, and 
Stop itself. Hence, it required very little user attention. 

It was desired that the GC operation should be as user 
oriented as the integrator operation. To achieve this some 
rather stringent constraints were placed on the use of the 
GC. The GC column was to be designed such that the five 
Gomponents? “GOP CCOSF Op" Nz seandsHlOrecoul dabemnescivedminio 
distinct peaks with a total analysis time of less than ten 
minutes. Further, no bypass valves, or other devices 


requiring user attention, were to be included in the column 
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design. 

Satisfying these constraints was not a»strivial task, as 
an examination of the gas chromatography literature quickly 
revealed that it is extremely difficult to easily separate 
this mixture into its individual components. By a trial and 
error procedure it was determined that the parallel-series 
column arrangement shown in Figure 6.7 would satisfy the 
stated constraints. From this figure it is seen that a 
Sample entering the column arrangement was equally split 
between the two parallel branches. The split ratio could be 
varied by adjusting the two needle valves located inside the 
GC oven, but it was found that an equal split gave the best 
results. 

Samples which passed through the Porapak T - Spherocarb 
branch of the column arrangement were resolved into three 
very distinct peaks, namely, a composite peak composed of 
Qe iemoanc1CO, a Gystinct; CO; ‘peak, anda ‘slightly taiding 
H,O0 peak. A typical chromatogram for a sample passing 
through this side of the arrangement is shown in Figure 
Coca. 

Samples which passed through the 5A molecular sieve 
branch of the column were also resolved into three distinct 
peaks, namely, an O, peak, a N, peak, and a CO peak. A 
typical chromatogram for a sample passing through the 5A 
molecular sieve column is shown in Figure 6.8(b). No peaks 
fOr H,Ovor CO, emergemtrom this cofumny as" ubese 4wo 


components are very strongly adsorbed by 5A molecular sieves 


z ‘ ; 4 tia ; u : if F : 7 
4) 4 ee ; 7 
- ui ya 
ear hy 
Ati ¥ ¢ 
’ ie 
iy ae v1 


ae 204), Lae ae ey LA MTIR ASS anne Bee 

\ grunge] log eigesaeA cS Fara to otaambunite 
| aripey fe 1G Yi one oe EE gud? selaween, 
Son, Laide a eras ab ceutily prey | af ore oR neg” 
| peti mete) ce 7) GSS toa 


so ire ee B22 egrend™ 2 wwone J4ceqepnatte nis LO Zz 


. 2 > wa hese @ ; .coinéetseceos besarte’ 


pe ek ("> he Fhe. ie ons on sauce 4iqmee 7 
oi | oagct Col! sre ses git ~-suted a 
i. . of Bs . 5 ats oo b@@gets ve beraee 7 

a a | om usa ‘a4 \ eats 40903. 880) 28 Fy pees 38 
- atlunes ; 

fy | 3 Kort stty eeblewtt Qed aor ee aang 

| ‘ - 
“ay SHeiteneg hte mine O4 oti is eb _ 


J 4a 8 | Toile 5 aw ay EF ie | - 225 yey ia 
i ae "eli } of . 2 ate # ewe ; _ 
cri’ Se (j a seer Pye) »BD< a | . a & - AGS « oe _ 


_ 
ember | e iting wae Y 30 ohiw cits cuss cf 
-; ‘+ 
7 Ao) <a ; 
a ee sane: Gayrasy dzidy vey a 


- 


oa. at 


116 


JuSwebuerry uwnto) yderbojyeuworys sey 7°9 aanbtg 


(WW OL'O) (ww O2'0) 
qied01aydS 1 yedeiog 
"00 ae 
| 
ee S inner 
O9/°N/°O : 
x BAJCA 
fe Buijdwes 
‘JaJOwelg |PUIWON = S85 95 
410}99q W ZEO0'O sue SuWNIOD II'7 
OL OL 
M8GE sinjessdwa! usao 
<O 
°N 
OO 


(O°H “°N 0 ‘209 ‘O9) 
a|dwes 


(W O'€) sanais (W 9°) 
JE/NDBIOW YS eqn, Aidwy 


d cd 
a 6. % 
PT - . 
ar 
a er ee y' € 
- ( } -s m4 ; 
ca if 


Le 
ese igo’). 
fg PYCISO BL 


—8UC 
az 
338 


- 


Recorder 
Response 


Recorder 
Response 


< O,/N,/CO 
Composite 


COs 


7 2 3 4 5 6 
Retention Time (min) 


] 2 3 4 5 6 
Retention Time (min) 


|< | —_0,/N,/CO 


Recorder 
Response 


Figure 6.8 


Composite 


pee Or 


] 2 3 4 5 6 
Retention Time (min) 


Typical GC Chromatograms 

(a) Porapak T — Spherocarb Columns 

(b) 5A Molecular Sieve Column 

(c) Total Column Arrangement as in Figure 6.7 
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at 353 K (the oven temperature). These components were 
purged from the column on a regular basis by raising the 
column temperature to 450 K for several hours with helium 
flowing through®’the columns. A” comparison” of “the (a) and 
(b) parts of Figure 6.8 shows that some overlapping of the 
peaks would occur if the exit streams from the two columns 
were mixed before the TC detector. To eliminate this over- 
lapping of peaks, a carefully selected length of empty 
tubing was placed in front of the 5A molecular sieve column. 
This empty column provided a time delay of approximately 30 
seconds which enabled the three peaks from the 5A molecular 
Sieve column to emerge in the 'window' between the CO, and 
the H,0 peaks from the other branch of the arrangement. A 
typical chromatogram of a sample passing through the total 
column arrangement is shown in Figure 6.8(c). 

Once the column arrangement was finalized, the GC was 
calibrated using various mixtures prepared by the previously 
calibrated flow metering and mixing equipment. These cali- 
brations are included in Appendix D along with a procedure 
for determining the overall sample composition from a 
knowlédge of the component response factors and the percen- 


tage area under each curve in Figure 6.8(c). 


6.6.2 Infrared Spectrophotometer 
The investigation of the dynamic behavior of the 
concentrations of CO, and CO from the reactor was performed 


through the use of a Pye Unicam model SP3-200 (Ser. #291036) 
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ratio recording infrared spectrophotometer. A Pye Unicam 
model SP3-050 (Ser. #299184) data processor was used in 
conjunction with the IR to perform specialized operations 
such as accurate placement of the monochrometer for constant 
wavenumber scanning and conversion of transmittance to 
absorbance. 

In order to continuously analyze a sample stream, it 
was necesSary to construct a continuous flow sample cell for 
Chewmimwetiwo celica were constructed,. tne first. of which was 
Found © EO be unsativstactoryedue to it having aevery low 
Signal to noise ratio. This was caused by a short cell 
pathlength and low IR transmission through the cell due to 
incorrect cell geometry. These factors necessitated the 
Gonsuructionm, Ofwas second sample .cell, shown tin Figure 6.9, 
which had a pathlength twice that of the first cell, and was 
constructed such that it conformed almost exactly to the IR 
beam geometry. These improvements resulted in satisfactory 
performance of the IR as the signal to noise ratio was 
increased by approximately an order of magnitude over that 
Ghucinable with whew finer cells tin theraonstruction. ot 
these cells, stainless steel was used throughout, with 
SoOdrTuMmchilorade windows in the) first cell, and) calcmum 
fluoride windows in the second cell. 

A further large increase in the signal to noise ratio 
was achieved by purging the IR sample chamber with air from 
which the atmospheric CO, and H,O had been removed. These 


gases were removed by passing low pressure air through a bed 
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Infrared Spectrophotometer Sample Cell 


Figure 6.9 
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of Ascarite (A.H. Thomas Co.) to remove the CO., and a bed 
of indicating drierite to remove the atmospheric H,0. 

Calibration curves for the IR, shown in Figure 6.10, 
were made by feeding the sample cell with gas mixtures which 
were prepared using the feed metering and mixing equipment. 
From these curves it 1S seen that the IR is approximately 
ten times as sensitive to CO, as it 1s to CO. In preparing 
these curves the CO absorbance was always meaSured at 2172 
em ', and the CO, absorbance was always measured at 2355 
cm '. At these wavenumbers CO and CO, strongly absorb 
Mbrarea radiation without interfering with eachyotherms It 
is seen that the calibration curves are not linear and thus 
do not conform to Beer's Law predictions. 

In calibrating the IR, aS in normal operation, the 
temperature of the gases entering the sample cell was always 
very close to the ambient temperature in the laboratory. 
This was achieved by passing the gases through the inner 
tube of a small double pipe heat exchanger which had cooling 
water at ambient temperature flowing between the two tubes. 
This heat exchanger was primarily built as a precautionary 
measure to prevent hot gases from the auxiliary reactor 


(500 K) from entering the IR sample cell. 
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7. OSCILLATORY BEHAVIOR IN THE OXIDATION OF CARBON MONOXIDE 
Once the equipment described in the preceding chapter 
was constructed, it was decided to investigate the existence 
of oscillatory behavior in the oxidation of carbon monoxide 
over supported metal catalysts. It was expected that this 
reaction could display oscillatory behavior due to the 
previous observations reported in the literature. These 
previously reported observations are summarized ina litera- 
ture review, and following the literature review the 
remainder of the chapter is devoted to a description of 
experimental observations of oscillatory behavior for CO 
Oxidation over supported metal catalysts uSing the equipment 


described in the preceeding chapter. 


7.1 Literature Review 

Presented in this section is a short review of the 
published literature concerning experimental observations of 
oscillatory behavior in the oxidation of carbon monoxide. 
The review 1S restricted to this rather narrow area as any 
comprehensive treatment of the various investigations 
concerning the catalytic oxidation of carbon monoxide would 
fill several volumes, with another volume necessary to 
reconcile the many contradictory observations which have 
been made. In the following literature review, the 
references have been cited in approximately their chronolog- 


ical order of publication so that some appreciation may be 
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gained of the sequencial developments in the field. 

Rewappears that peusen et al. [7.1 7.2] were the: first 
to observe oscillations in the catalytic oxidation of carbon 
monoxide. Their observations were made using a single 
pellet of 0.3% Pt/Al,0, suspended in a differential reactor 
(homrecycle) “Oscillataons in ’the carbon) dioxide concentra- 
tion were observed in the temperature range 453 - 523 K, and 
typically these oscillations were accompanied by small ( ~ 
2 K) temperature oscillations. They speculated that several 
states of chemisorbed carbon monoxide participated in the 
reaction. This hypothesis was further refined by Hugo and 
Jakubith [7.3] who observed carbon dioxide oscillations in 
the temperature range 384 - 423 K when using a platinum mesh 
catalyst in a recycle reactor. They proposed that the 
oscillations were the result of a shift between bridged and 
linear forms of chemisorbed carbon monoxide. 

In the two preceeding works, observations of oscilla- 
tions were made by analyzing gas phase compositions, thus 
oscillations on the catalyst surface could only be inferred. 
The existence of surface oscillations was first estabilished 
by Dauchot and Van Cakenberghe [7.4] who showed that concen- 
tration and temperature oscillations could occur on catalyst 
surfaces. This was done through the use of resistivity 
measurements for platinum wires, and work function measure- 
ments for thin platinum films. Their conjecture was that 
the oscillations were caused by differences in the rates of 


adsorption of oxygen and carbon monoxide, with large surface 
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thermal effects resulting from the adsorption and reaction 
processes. 

Bach of the preceeding works had in common the use of a 
platinum catalyst and the rejection of the classical trans- 
port phenomena explanation for the oscillatory behavior. 
Bothgof.theseslinkspwere broken in, work by Eckert .et al. 
[Jeolein whichsconcentration and temperature oscillations 
(~3 K) were reported for carbon monoxide oxidation in a 
Becyc lesneac COmeIOeNhECh eat CUOsAlZ On with) jcatalyst was 
used. The recycle rate and the reaction temperature (430 K) 
were deliberately set such that transport limitations were 
encountered, and this was used as the explanation for the 
oscillatory behavior. 

However, the earlier view that oscillatory behavior 
could occur in the absence of transport limitations was 
reinforced by McCarthy et al. [7.6] who showed that oscilla- 
tory behavior could occur in a gradientless, spinning 
ba Stet wiGOme VNONS Starred tank wreactorms sAm0eCS5 tert Ado. 
catalyst was used and the oscillations were observed at 
reaction rates near the maximum of the 'abnormal' reaction 
rate curve for carbon monoxide oxidation. It was believed 
that these oscillations were due to nonunique surface reac- 
tion ratescontrol salt subsequent wonkebyweVenghese se laa in. 
[7.7] using the same rector and catalyst in the temperature 
range 37.0 4 430.K, it wwas found sthat byvo wuelitativery 
different types of oscillations could occur depending on 


whether ‘pure’ or ‘impure' oxygen was used. When oxygen 
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with trace amounts of hydrocarbon impurities was used, large 
multipeak and chaotic relaxation oscillations were observed. 
However, if very pure oxygen was used only small symmetrical 
oscillations occurred. The existence of the relaxation 
oscillations was attributed to periodic build-up and 
elimination of the hydrocarbon impurities on the catalyst 
surface. No explanation was attempted concerning the small 
Symmetrical oscillations, but the existence of localized 
temperature excursions could not be ruled out as a possible 
cause. 

In parallel with the above work, oscillatory behavior 
resulting from the simultaneous oxidation of hydrocarbons 
(1-butene) and carbon monoxide was explicitly studied by 
Cublapeand@kenney (7.160) Usingmal 0s5%ePtyYAl, Osecatalyst in a 
gradientless recycle reactor. They failed to observe any 
oscillatory behavior when only carbon monoxide was oxidized, 
but suStained oscillations and multipeak behavior were 
encountered in the temperature range 423 - 438 K when the 
Simultaneous oxidation of 1-butene and carbon monoxide was 
Carried out. They believed that the oscillations were the 
Fesulb@ol al périodice#bui ld-uprand burneoiteci athe —purere 
onbthemcatalystrsuntacer 

The belief that multipeak and chaotic phenomena were 
the result of the interaction of two reactions (hydrocarbon 
and CO oxidations) was challenged by Sheintuch and Schmitz 
[esi eePlichtawancdescnimtz [7.100 and@Sheimeuch [(7.71)° in 


investigations of carbon monoxide oxidation on a platinum 
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foil in a gradientless stirred tank reactor. They found 
that single and multipeak behavior, multiple limit cycles, 
and chaos could occur in the temperature range 423 - 544 kK. 
It was their belief that intrinsic chaotic behavior exists 
for this reaction. However, this assertion must be 
tempered with the observation that the activity of the 
platinum was not very aeeute Wake. 

The existence of complex dynamic phenomena has also 
been reported by Rathousky et al. [7.12] in an investigation 
ObecalLbonemonoxider oxidation on avPt/Al.O, ‘catalyst in a 
recycle reactor over the temperature range 433 - 473 K. 
They present experimental results which indicate that long 
periods of time (>9 hours) are necessary to determine the 
character of some types of complex dynamic behavior. 

In summarizing the experimental work which has just 
been described about all that can be said is that oscilla- 
tions do exist, and these oscillations may be simple or 
complex in nature. No concensus has yet been reached 
concerning the underlying mechanism which causes the oscil- 
lations, and in fact very little work aimed at elucidating 
this mechanism has been performed since most of the reported 
experimental work has been of a descriptive (as opposed to 
mechanistic) nature. However, despite the abundance of 
descriptive studies, a thorough experimental examination of 
the effects of various operating conditions on oscillatory 


behavior has yet to appear in the literature. 
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In the following sections of this chapter, a detailed 
investigation of the effects of various reactor operating 
conditions on oscillatory behavior is presented. Whenever 
possible, comments are made regarding the underlying 


mechanism responsible for the oscillatory behavior. 


7.2 Preliminary Experimental Work 

At the time that the experimental work was started it 
was not possible to state a priori whether or not a 
Particular catalyst could display oscillatory behavior. 
Hence, a certain amount of preliminary work with several 
catalysts was necessary in order to determine both the reac- 
tion characteristics of the catalysts as well as the 
operating characteristics of the apparatus. 

The first catalyst which was examined was a supported 
platinum on alumina catalyst obtained from Englehard Indus- 
tries Limited. This catalyst, which contained 0.3 weight 
percent platinum, “was obtained in the slorm. of dcylindrical 
alumina pellets which were one-eighth of an inch in diameter 
and length. Several beds of this catalyst were placed in 
the reactor and used for preliminary investigations into the 
Macoure, of the catalytic oxidation Olwcarbon monos Ochmmenn 
each case the catalyst pellets were pretreated in flowing 
Grywair at 625K ineastube Lurmace fcrgapproximave ly eewel ve 
hours before being placed in the reactor. Also, the weight 


of the bed placed in the reactor was kept constant at twenty 
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grams. If less than twenty grams of catalyst pellets were 
used, then the remainder of the bed was composed of inert 
alumina cylinders which were similar in size and shape to 
EnexecatalystipelietsemeUntortunately, «thesguse sof “this cata- 
lyst did not lead to any observations of oscillatory behav- 
10r. Instead it was found that this catalyst displayed 
other interesting types of behavior. In particular, at low 
temperatures (363 - 373 K) this catalyst displayed multiple 
Steady-state behavior in that for a given set of operating 
conditions it was found that two stable steady-states always 
existed. These steady-states represented the extremes which 
wenempOossible;t:.erutotalaconversion,; or virtually no 
conversion of the carbon monoxide. The high conversion 
Steady-State could be reached by initially purging the 
reactor with pure oxygen before introducing any carbon 
monoxide in the feed stream, and conversely the low conver- 
Sion steady-state could be reached by initially purging the 
eae erOrewht hipuneskearbon monoxide (or sang dA %"CO, in N, 
mixture) before introducing any oxygen in the feed. Also, 
it was sometimes possible to drive the system from a high 
conversion steady-state to a low conversion steady-state by 
making an abrupt increase in the carbon monoxide concentra- 
tion in the feed. However, the reverse was not true as a 
low conversion steady-state could not be driven to a high 
conversion steady-state by making an abrupt decrease in the 


feed stream carbon monoxide concentration. 
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It was with regard to these transitions from high to 
low conversion steady-states that a further interesting 
property of this catalyst was observed. It was found that 
the previously mentioned transition from a high to a low 
conversion steady-state did not occur if the carbon 
monoxide concentration was increased incrementally rather 
Ehammeaueupuly jean ntact oi hethercarbon monoxide sconcentra> 
tion waS increased incrementally it was possible to feed 
very large amounts of carbon monoxide to the reactor and 
Still stay at the high conversion steady-state. However, as 
the molar flow of carbon monoxide became large the reactor 
ceased to behave isothermally, as small (2 - 5 K) tempera- 
ture gradients developed over the length of the catalyst 
bed. 

This type of behavior (total conversion of large 
amounts of CO with a temperature gradient along the bed) 
persisted even when only very small amounts of catalyst 
(less than one gram) were placed in the reactor. It was 
believed that this behavior was the result of catalyst 
surface temperatures which were much higher than the bulk 
temperature. If a finite amount of time was required for 
these surface temperatures to adjust to feed concentration 
changes, then this hypothesis also explains why the very 
high reaction rates could only be reached by incrementally 
increasing the feed carbon monoxide concentration. 

After investigating the behavior of this catalyst at 


many different operating conditions it became apparent that 
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oscillatory behavior was not going to be displayed, and that 
it would be necessary to use a different catalyst. 

The second catalyst which was used in these preliminary 
screening tests was a General Motors converter catalyst 
(type #78925, package #6498995). This catalyst, which 
contained approximately 0.04% platinum and 0.04% palladium 
Supported on alumina, was obtained in the form of spherical 
pellets approximately one-eight of an inch in diameter. 

This catalyst was primarily chosen because it had a lower 
metal loading than the first catalyst. “By using a catalyst 
with a low metal loading it was hoped that isothermal 
reactor operation could be maintained, and that some of the 
problems encountered with the first catalyst could be 
avoided. 

After investigating the behavior of this Pt-Pd catalyst 
at many sets of operating conditions, it was found that 
oscillatory behavior ycould occur, and@shown in Figure 7.1 
are some typical oscillations of the carbon dioxide concen- 
Brea cloOnmrase CCOLdecaaby using the ile lat. 25 S0ucnet ss SlOwnien 
Figure 7.1(a) is a very long period oscillation which 
occurred when five grams of catalyst (diluted to twenty 
grams) were placed in the reactor. This oscillation was 
very regular as it was reproducible over a period of several 
hours. It was found that the period and amplitude were 
guite sensitive to the concentration of CO, ingthe feed, .and 
iti wasepossible tolobtain oscillations;with periods in 


excess of one hour. 
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(a) T = 413 
(b) T = 413 
(c) T = 403 
Cd) oT = 403 
(e) T = 403 
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Since studying the parametric sensitivity of an oscil- 
latory phenomena with a period of an hour would be extremely 
time consuming, the operating conditions and the amount of 
catalyst in the reactor were varied in an attempt to obtain 
oscillations with shorter periods. Shown in Figures 7.1(b) 
and (c) are oscillations which occurred when twenty grams of 
the catalyst was placed in the reactor (no dilution). These 
oscillations were recorded using the first IR sample cell 
which had a low signal to noise ratio as mentioned in the 
previous chapter. The increase in the signal to noise ratio 
which was achieved through the design of the second IR 
sample cell is shown in Figure 7.1(d). 

It 1S seen that what appears to be a single peak oscil- 
lation in Figure 7.1(c) could in fact have been a two peak 
pew cycléevosc¥llation as shown. in-Figure 7.1(d}); as the 
small peaks which are readily observable using the second IR 
sample cell would have been completely obscured by the noisy 
Signal which was obtained from the first IR cell. 

Shown in Figure 7.1(e) is another interesting phenomena 
which was observed during the preliminary work with the 
twenty gram bed of Pt-Pd catalyst. While an oscillation was 
definitely occurring, it is seen that the behavior was not 
periodic as there was no regular cycle which repeated 
PeSeLi.. This typegoiiebehavior was, foundutc occur governa 
fairly large range of operating conditions. 

Once it had been established that oscillatory behavior 


déhwnieely soccumreaguwhen®this Pt-Pdycatalyery was used, it 
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was decided to determine if oscillatory behavior could be 
observed if a different catalyst was used. To this end a 
number of supported platinum catalysts were prepared. These 
catalysts contained 0.05 weight percent platinum supported 
on one-eighth inch alumina spheres. Two different batches 
of this catalyst were prepared, and although both batches 
displayed oscillatory behavior, they did so under different 
reaction conditions due to large differences in activity 
between the two batches. 

Tiemrirst. batch orm. 05% PtyvAIs Oe catalyst, was not 
homogeneous as large dark spots, indicating high metal load- 
ings, were visually apparent on many of the catalyst 
particles. While this catalyst displayed oscillatory behav- 
ior, as shown in Figure 7.2(a), the nonhomogeneous nature of 
the catalyst resulted in a very erratic conversion of CO, 
with a resulting high noise level in the IR signal. Also, 
it was suspected that local hot-spots existed on the cata- 
Iyst=particies during reaction conditions due to the high 
reaction rates in the regions of high metal loadings. 

Some of these problems were alleviated when a second 
batchoft 10.05% Pt/Al.O, catalyst was prepared using da rotany 
Grier instead of manual stirring during the evaporation step 
in the catalyst preparation. These catalyst particles were 
visually very homogeneous, and it was found that the noise 
level in the IR signal was much reduced when this second 
batch of catalyst was used instead of the first batch. 


Also, the second batch of catalyst had a much higher 
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acervrtyethan@either the first batch; or the Pt-Pdi converter 
catalyst, as it displayed oscillatory behavior at very low 
temperatures at which the other catalysts were virtually 
piactune..) SHOWN In Figures 7.2(bp) to (d) are oscillations 
which were observed when twenty grams of the second batch of 
0.05% Pt/Al,O, catalyst were used in the reactor at 318 K, 
S2oeeacidasoo Kk, crespecuively., Wt was found that the rela-— 
tively high noise level apparent in Figure 7.2(d) could be 
reduced by increasing the volumetric flow rate, as shown in 
Buigures 72. e):, 

Theprosci lations shown. in Figure 7.2 are typical of 
Enosemobservea for the two batches of°0:05% PE/AL.O, cata- 
lyst, aS it was rare that a periodic single peak per cycle 
oscillation would develop. Rather, a very irregular oscil- 
facvon woulcsusually occur, and typically this oscilbation 
would not continue for more than a few hours. Usually 
within a few hours (sometimes a few minutes) of the start of 
the oscillation, a large CO maximum (CO, minimun) would 
occur which would drive the system to a low conversion 
Steady-state. 

Thus, in the preliminary work it has been shown that 
oscillatory behavior can occur for several different cata- 
lysts. Of the catalysts used, the Po sPapconververicacaiyca 
gave the 'best' results in the sense that besides chaotic 
behavior, a periodic oscillation with a low noise level 
could be observed. The existence of a periodic, low noise, 


oscillation was desired so that the effects of changes in 
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operating conditions on the period and amplitude (or the 


existence) of the oscillatory behavior could be examined. 


7.3 Effects of Operating Conditions 

In this section are presented the results of a detailed 
eXamination of the effects of various operating conditions 
on the existence and behavior of oscillations observed for 
the oxidation of carbon monoxide carried out on the Pt-Pd 
converter catalyst described in the previous section. All 
of the observations presented in this section were obtained 
MSsimgmewenty tans OL une, GM convertor catalyst, and in all 
cases the carbon monoxide was obtained from a specially 
prepared mixture of 11% carbon monoxide in nitrogen. 

The oxygen which was used was normally Linde very dry 
grade, but this was replaced by Matheson ultra high purity 
oxygen on a number of occasions, and no noticable effect on 
the oscillatory behavior was observed. Also, it was deter- 
mined that the reactor, recycle pump and associated tubing 
were catalytically inert by observing that no carbon dioxide 
was produced when carbon monoxide - oxygen mixtures were 
Bonet ed to the reactor system at 425 K with the catalyst 
bed removed from the reactor. Unless otherwise stated, the 
reactor pressure was always set at one psig so that the 
absolute reactor pressure was very close to one standard 


atmosphere (compensation for Edmonton's elevation). 
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ae COn COs andy Combined: OSci l lations 

Shown, ineFigures 793 a) and (b) are typical oscilla- 
PLOncmOLeCO, candyCO respectively, at 393°K. ‘These cscivila- 
tions were both self-starting and self-sustaining in that no 
external perturbations were necessary in order to have the 
oscillations start, and once started the oscillations 
continued indefinitely. The oscillations shown in Figures 
7.3(a) and (b) were not recorded simultaneously since the 
IR spectrophotometer can only monitor one frequency at a 
time. The slight mismatch in the oscillations is due to the 
Sequential nature of the recordings. 

Mnesoserllatvon ine carbon dioxide concentration) shown 
in Figure 7.3(a) can be broken down into three regions which 
continually repeat themselves, namely, a slow increase to a 
maximum, a more rapid decrease to a minimum, anda fast 
(virtually instantaneous) increase to a local maximum with a 
Small amount of overshoot. These three regions had counter- 
parts in the carbon monoxide oscillations shown in Figure 
fests) which were: a region of vartually “zero CO concentra. 
tion, a slow rise to a maximum, and a rapid decrease to zero 
concentration, respectively. 

Shown in Figure 7.3(c) are the oscillations which 
resulted from passing the effluent from the main reactor 
through the auxiliary reactor such thatralis thescanbon 
monoxide was totally converted to carbon dioxide. The 
approximately straight tracing in the rightmost part of 


Figure 7.3(c) is the result of bypassing the feed stream 
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Figure 7.3 Simple Oscillatory Behavior 
(a) Periodic Carbon Dioxide Oscillations 
(b) Corresponding Carbon Monoxide Oscillations 
(c) Effluent From Auxiliary Reactor 
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around the reactor such that it passed directly through the 
auxiliary reactor, thus this tracing represents the 100% 
conversion line. From a comparison of Figures 7.3(a) to (c) 
it 1s seen that there is no phase lead or lag between the CO 
ana cOs- concentration oscimratvons. as the oscillations in 
Figure 7.3(c) appear to be a direct combination of those in 
Figures 7.3(a) and (b) once the differences in IR 
sensitivity to CO and CO, are taken into account. If there 
was any phase lead or lag between the CO and CO, oscilla- 
taons, lt would show up imo the form of local maxima, minima, 
eno ,OreextraspCintus sof sinmlection wins hagures7 23 (cyer This 
absence of phase differences between the CO and CO, oscilla- 
tions can be used aS an indication that the oscillations are 
occurring in the absence of mass transport limitations. 

A detailed examination of the dynamic behavior of CO,, 
GG, -anaecO, sphus COsshown wn Figure+?.3s yields valuable 
insight into the phenomena occurring during oscillations. 
MhemesGiiltctioneat aboutei? to 77 imuneinGckiguret sy ssewill ibe 
used as an illustration. (It should be kept in mind that 
the slight mismatch in oscillations between parts (a), (b) 
and (c) of the figure are due to the previously mentioned 
sequential recording of the outputs.) Observations for the 
abovetoseitlation, which*are?’similat to those foutouner 
oscillations, are summarized in Table 7.1. 

The observations in Table 7.1 are explained by the 
following sequence of events: In Region 1, CO accumulates on 


the catalyst surface and a gradual increasé in CO, “rate of 
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Taptes yes Descriptionsot Oscillation Shown ain Figure’ 7.3 


Region Time Behav LOrmot, BulkeGOn. COuand 


COsplus CC> Concentrations 


“ecsadvalslincrease in<COm conc: 
Coj)eeenou COs in product 
(CUREconDinedmcosp lis ICO. Bin forodtc t 
less than CO in feed 


2 (A OmcOe OO sue) mec an hy svaprorincreaséeé ineCO, conc. 
(os farmiverapia increase: in CO <cone). 
Ci ecombineducO plus sCO, in product 
exceeds amount of CO in feed 


3 Lomo Guia Cay veryurapiamincrease in CO. conc. 
(bb) very. raped decrease in CO conc. 
COBZ EES 


(cormecombinedecOnplis™CO7, inv product 
decreases from greater than to 
less than amount of CO in feed 


formation occurs; in Region 2; desorption of CO ana 
decreased rates of CO, formation take place; and in 

Region 3, rapid accumulation of CO on the surface of the 
Catalyst and “rapid Sincrease in CO, production ‘occurs: 
Significant oscillation in the temperature of the metal 
crystallites appear to be the only reasonable explanation 
for these phenomena. The increase in CO accumulation and 
CG, formatiomein Region “results! inwan anegrease! ini ‘surface 
(or metal crystallite) temperature. This increase has to be 
very rapid at the end of Region 1, which results in the CO 
déesorpiiom in Region 2.8 ‘The! cooling “ebLectsi/of «CO mesorp-— 
tion tend to balance the heating effects of CO reaction 


until the surface is completely depleted of CO. Once the 
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surface is depleted, the temperature drops rapidly with 
Subsequent: rapid aesorptmon wf CO andi increased CO. Sproduc- 
tion in Region 3. The mechanism by which communication 
between different metal crystallites in the catalyst pellets 
occurs is not clear. Changes in bulk CO concentration 
and/or slight variations in bulk phase temperatures may be 


responsible. 


7.3.2 Recycle Ratio 

Shown in Figure 7.4 are the effects of recycle ratio on 
aucanponsdroxide oscillation. It Wis seen in Figure 7.4a) 
that if the reactor is operated as an integral bed reactor 
inomnecyclejy that oscillatory behavior canioccurw, Upon 
close examination the oscillations from the integral bed 
reactor are seen to be not simple in nature, but rather each 
cycle is composed of two distinct peaks. This observation 
is reinforced in Figure 7.4(b) where it is seen that the 
positive feedback effect which results from increasing the 
recycle ratio from 0:1 to 3:1 has resulted in the doublet 
nature of the oscillation becoming emphasized. However, as 
the recycle ratio is further increased in Figures 7. 44c) to 
Chyeite rs seen that a new form of OSGi ldatronmcn miaaamto 
that shown in Figure 7.3 emerges. The oscillation shown in 
Figure 7.4(f£) did not change upon further increases in the 
recycle ratio and hence this oscillation represents that 
which occurs in a well mixed reactor in the absence of 


transport limitations. 
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It 1S believed that the oscillations shown in Figure 
7.4(a) were caused by transport resistances between the bulk 
phase and the catalyst surface, whereas the oscillations in 
Figure 7.4(f) occurred in the absence of transport resi- 
Stances, and thus were possibly related to the underlying 
mechanism of the reaction. Since it was desired to examine 
oscillations in the absence of heat cee TARE transfer limi- 
tations, the recycle pump was always run at 1800 RPM which 


corresponded to a recycle ratio of between 100:1 and 200:1. 


7.3.3 Bulk Temperature 

Shown in Figure 7.5 are the effects of bulk temperature 
on the oxidation of carbon monoxide. It is seen that as the 
bulk temperature was increased, the system moved from an 
intermediate conversion steady-state shown in Figure 7.5(a), 
to a periodic single peak per cycle oscillation as shown in 
Figure 7.5(b). Upon further increases in the bulk tempera- 
ture it is seen in Figure 7.5(c) that the frequency and 
amplitude of the oscillations decreased, and that the oscil- 
lations were no longer periodic, as shown in Figures 7.5(d) 
and (e’). As the bulk temperature was further increased, the 
frequency and amplitude of the oscillations became very 
small as shown in Figure 7.5(f£), and finally the system 
passed to a high conversion steady-state as seen in Figure 
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7.3.4 Concentration of CO in the Feed 

Shown in Figure 7.6 are the effects of varying the 
concentration of carbon monoxide in the feed to the reactor. 
It 1s seen that decreases in the feed concentration of CO 
resulted in the system behaving quite similarily to that 
Shown in Figure 7.5 for increases in the bulk temperature. 
In particular, as the feed CO was decreased, the system 
passed from an intermediate conversion steady-state to a 
periodic single peak per cycle oscillation. Upon further 
decreases in the feed CO concentration, it is seen that the 
frequency and amplitude of the oscillations decreased, and 
ultimately the system passed to a high conversion steady- 
State at a very low concentration of CO in the feed stream. 

It 1S interesting to note that there was also a quali- 
tative change in the shape of the oscillations as the feed 
CO concentration was lowered (or the bulk temperature 
increased), as the slow increase to a maximum after the 
rapid increase in CO, concentration which is clearly shown 
in Figure 7.6(b) is virtually eliminated in Figures 7.6(d) 
Atlmeues) o)) LE 1s belseveduthat this) Slow imicreases anetnemce. 
concentration is the result of a surface capacitance effect. 
Hence, it would be expected that this effect would be more 
prominant at high concentrations of €QJin thesfieed, as this 
should correspond to high surface concentrations of carbon 
monoxide. Similarily, low bulk temperatures should imply 


large CO surface concentrations. 
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7.3.5 Volumetric Feed Rate - CO/O, Constant 

Shown in Figure 7.7 are the effects of changes in the 
volumetric feed rate on the behavior of the system. In 
examining these effects the bulk temperature and the ratio 
of CO to O2. were held constant. It is seen that decreases 
in the volumetric feed rate resulted in qualitatively ident- 
ical behavior to that shown in the preceeding section for 
decreases in the feed CO concentration. As shown in Figure 
7.7, the system passed from an intermediate conversion 
steady-state to a large amplitude periodic oscillation when 
the volumetric feed rate was decreased. Upon further 
decreases in the volumetric feed rate, the frequency and 
amplitude of the oscillations decreased, and it is believed 
that ultimately a high conversion steady-state would have 
been reached. This waS not attempted as it would have 
necesSitated going to very low feed rates with a resulting 
large increase in the noise in the IR signal due to poor 


mixing in the IR sample cell. 


7.3.6 Volumetric Feed Rate - Molar CO Flow Constant 

In Figure 7.8 are shown the effects of variations in 
volumetric feed rate to the reactor with the molar flow of 
carbon monoxide held constant. It 1S Seen that as the volu- 
metric feed rate was increased that the system passed from 
an intermediate conversion. steady-state to a periodic single 
Deak per cyClemosccultacion, » HOweVCigmocethesteed raterwas 


further increased, a comparison of Figures 7.8(b) and (c) 
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shows that the frequency and amplitude of the oscillations 
were virtually unchanged. A further increase in the volu- 
metric feed rate resulted in a decrease in the frequency and 
an increase in the amplitude (using the conversion scale) of 
the oscillations as seen in Figure 7.8(d). Further 
increases in the volumetric feed rate were not made, as 
further dilution of the carbon monoxide and carbon dioxide 


resulted in very weak, noisy, signals from the IR. 


7.3.7 Nitrogen and Helium Diluents 

Shown jin Figure 7/9vare the effects on the carbon 
monoxide oscillations due to the partial replacement of 
oxygen with nitrogen and helium. The oscillations 
preceeding point A in Figure 7.9 were typical of those which 
occurred when the reactor feed was composed of carbon 
monoxide and oxygen (with a small amount of nitrogen from 
the CO/N, mixture). At point A, nitrogen was used to 
replace approximately 25% of the oxygen in the feed and it 
is seen that this had a definite effect on the oscillations. 
A comparison of the oscillations before and after the intro- 
duction of the nitrogen shows that the widths of the oscil- 
latory maxima were increased due to an elongation of the 
regions between the maxima and the abrupt transitions to 
100% conversion. 

At point «B omn *Fagure (69 the ver cde wivone weomitcro- 
duced at point A, was replaced by an equal flow of helium. 


It is seen that replacing the nitrogen with helium had only 
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a minor effect on the carbon monoxide oscillations, as the 
small differences before and after point B can be attributed 
to inaccuracies in the flow system by which nitrogen and 
helium were used to replace oxygen. Thus, it appears that 
the changes in the oscillations which result from replacing 
oxygen with nitrogen or helium are solely due to a dilution 
effect, as equal dilutions with nitrogen and helium had 
equal effects on the carbon monoxide oscillations. 

These observations seem to rule out heat transfer from 
the metal crystallites to the gas phase being important to 
the existence of oscillatory behavior, as otherwise the 
large difference in thermal conductivity between nitrogen 
and helium would be expected to have a more pronounced 
effect than was observed. 

Re DOM tern Inehiaure 739 cas further 25g) or the oxygen 
feed was replaced with helium, and it is seen that this had 
a very dramatic effect on the oscillations, as the regions 
between the maxima and the abrupt transitions to 100% 
conversion became very elongated. It is worth noting that 
these were the only parts of the oscillations which are 
seriously affected by the replacement of oxygen with helium, 
as the abrupt transitions to 100% conversion, the flat 
region at 100% conversion, and the slow rise to a maximun 
were all only slightiy changed from those originally (0)— 40 
minutes in Figure 7.9) observed. This seems to indicate 
that the oxygen concentration is only of importance during 
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factors (surface temperature, CO concentration, etc.) are 


important during the other parts of the cycle. 


7.3.8 Reactor Pressure 

As described in the previous chapter, the reactor pres- 
Sure could be varied by changing the external set point 
pressure applied to the reactor pressure controller. Shown 
in Figure 7.10 are the effects of varying the reactor pres- 
Sure from 1 to 16 psig (1 to 2 atmospheres absolute). It is 
seen that the amplitude of the oscillations continually 
decreased as the pressure was increased, but that the 
frequency of the oscillations passed through a minimum. 

This type of effect in which the frequency passed 
through a minimum when an operating parameter was changed 
had not been previously encountered, and hence these obser- 
vations were carefully checked and they were found to be 
reproducible. As the pressure 1S increased, it appears that 
competing processes occur which affect the surface concen- 
trations of oxygen and carbon monoxide due to changes in the 
bulk phase concentrations. However, the manner in which 
these processes produce the effects shown in Figure 7.10 is 


not presently understood. 
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7.4 Qualitative Effects of Temperature 

In most of the preceeding sections, the reaction 
Penveralure was sheldsoonstant at 9393 Ke to; facilitate compar- 
lons between the effects of the various operating conditions 
which were varied. However, it has been found that inter- 
esting changes occurred in the oscillatory behavior when it 
was examined at temperatures other than 393 K. Some of 
these qualitative changes which occurred when the bulk 
temperature was varied are shown in Figure 7.11, where it is 
seen that two qualitatively different types of oscillations 
could occur depending on the reaction temperature. 

It is believed that different mechanisms are respons- 
ible for the different types of oscillations shown in Figure 
Uae) Musue ek); andud. Ce) (423K) 2° However, ithe 
mechanisms are obviously not mutually exclusive as shown by 
the coexistence of both types of oscillations in Figure 
7.11(d) (413 K). As was stated when the effect of the 
recycle ratio was examined, it is believed that the small 
fast oscillations are the result of transport resistances, 
while the large relaxation oscillations are probably related 
to the underlying mechanism of the reaction. 

Further evidence that transport resistances are signif- 
icant at 423 K can be inferred from the steady-state 
behavior of the reaction. At temperatunes Tess than) 213% 
no high conversion steady-state was ever observed which was 
sroniticantly+duéterent from +l00% *convercicivje).c. ane sports 


could ever be obtained on the reaction rate curve in the 
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region where the reaction rate increased with increasing 
Garoone monoxide concentration, as Ehe CO*concentration was 
never significantly different from zero when the system was 
at a high conversion steady-state. Thus at temperatures 
less than 413° K the reaction rate was not limited by any 
transport constraints, and the reaction rate controlling 
step was one of the kinetic steps in the reaction. 

However, at temperatures above 413 K high conversion 
Steady-states at less than 100% conversion could be 
obtained, as at these steady-states it was possible to 
measure finite amounts of carbon monoxide in the reactor 
effluent. Since these high conversion steady-states were 
normally obtained under conditions of high CO concentration 
in the feed stream, it was suspected that what was being 
measured was a transport resistance curve, and not a reac- 


tion rate curve. 


7.5 Chaos and Nitrogen Dilution at 423 K 

In the preceeding section it was shown that as the 
temperattiré-incréased the large* relaxation Oscillations 
disappeared and only the transport related oscillations 
remained. Since increased rate of reaction is the main 
result of increased temperature, it was suspected that the 
pelaxation oscillations could be observed at high tempera 
tures if the reaction rate was somehow decreased. This 


decrease in the reaction rate was most easily accomplished 
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by replacing oxygen with nitrogen in the reaction gas 
Meeciremate47enhKandatheselfects Of this sat, various 
nitrogen to oxygen ratios are shown in Figure 7.12. 

It 1S seen, as the nitrogen to oxygen ratio was 
increased, that the occurrence of relaxation oscillations 
could be initiated. Further increases in the nitrogen to 
oxygen ratio resulted in the frequency and amplitude of the 
relaxation oscillations increasing until they completely 
obscured the transport related oscillations. Also, it is 
seen that the oscillatory behavior was not regular and 
periodic, as in the preceeding sections, but rather was very 
irregular, and perhaps should be called chaotic. 

As shown in Chapter 5, when two nonlinear oscillators 
are operating simultaneously, interesting phenomena such as 
@haos,can occur. Thus, “the chaotic behavior which was 
observed for this catalytic system 1s attributed to the 
interaction of the mechanistic and the transport related 


SSCHisations,. 


7.6 Reproducibility of Oscillations 

Se RGEE in the catalyst activity can be a major problem 
when examining catalytic reactions, for if the activity 
changes then experimental observations are not reproducible 
and very little, information Gan be Obtained vabour sw hesinue 
mechanism of the weaction. “Fortunately, In this study.) it 


Wes sound thatethe catalyst activity was Very constant. 
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This 1s shown in Figure 7.13, where it is seen that the 
oscillatory behavior was very reproducible over a period of 
approximately 250 hours of use of the catalyst. 

inerigure 7.13), thevsoscillations shown in the (a) 
section were obtained at the beginning of the preceeding 
investigation of the effects of operating conditions on 
Gseillatory behavior, and the oscillations in the (b) and 
(c) sections were obtained in the middle (after ~150 hours 
of operation), and at the end (after ~250 hours of opera- 
tion) of this investigation, respectively. While these 
oscillations were very reproducible, it was found that 
changes in catalyst activity could occur under certain types 
of reaction conditions. | Fortunately these particular reac- 
tion conditions were not encountered during the investiga- 
tion of the effects of operating conditions, but rather were 
encountered while performing subsequent work at higher 
Shy Sea eRe pee 

These changes in activity are well illustrated in 
Figure 7.14 where the effect of severe reaction conditions 
at 473 K on the catalyst activity were investigated. It is 
seen that interesting transient phenomena Occurred abcer 
each increase in the CO feed concentration at points A, B 
and C, and it appears that a sustained oscillation was 
occurring between points C and D. At point D the reed CoO 
concentration was further increased, and as can be Seen, an 
abEupte transition vo 4a Low GOnVersion Steady-state occurred. 


After the system moved to the low conversion steady-state, 
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the feed CO was stopped, and the reactor was purged with 
oxygen. 

Purging with oxygen resulted in the carbon monoxide 
concentration in the reactor being decreased to the point 
where the oxidation reaction once again started to proceed 
with an appreciable rate. In fact, the rate of the reaction 
became so large that the CO concentration in the reactor 
decreased to a low level virtually instantaneously. At 
point E, the reactor feed conditions were reset to those 
which existed before point A, and changes in the reactor 
feed conditions at points F, G and H were made which were 
Similar to those previously made at points A, B and C, 
respectively. 

A comparison of the system behavior after the similar 
operating condition changes were made at points A and F, B 
andeGeeorn Cand H-showssthat osignifticant’ changes occurredsin 
the catalyst activity. It is apparent that the catalyst 
activity was much higher after point E, than it was before 
point D. It is believed that the increase in activity 
probably occurred due to some restructuring of the catalyst 
surface which occurred when the system was moved from the 
low to the high conversion steady-state by purging with 
oxygen. This restructuring was possibly caused by extremely 
high surface temperatures which occurred when the reaction 
rate became very large as the CO concentration in the 
ceactorswentefrom a high toua low level virtually 


instantaneously. 
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To verify that the catalyst activity had indeed 
changed, the reactor temperature was reduced to 393 K, and 
the feed conditions were set at values similar to those 
wivehwine sulted sin “sheroses lationseshown ‘insFigure wisl3ey It 
was found that oscillatory behavior no longer occurred as 
the system went to a total conversion steady-state due to 
Ene wincreased activity of sthe catalyst. Thus, the oscilla- 
tions shown in Figure 7.13 could no longer be reproduced 
after subjecting the catalyst to the severe reaction condi- 
ErOnSere tet 7S ok 

These observations of changes in the catalyst activity 
were of some concern, as they seemed to throw into doubt the 
reproducible nature of the oscillatory behavior. An attempt 
was made to determine why the catalyst behavior had origin- 
ally been so reproducible. From an examination of the 
records of the preliminary work with this catalyst, it was 
determined that the catalyst temperature had never been 
raised above 433 K under reaction conditions (625 K during 
pretreatment). Also it was determined that at 423 K the 
catalyst had several times been subjected to the previously 
mentioned severe reaction conditions, i.e. at 423 K the 
reaction several times underwent abrupt transitions from low 
to high conversion steady-states due to the reactor being 
purged with oxygen after it had filled with carbon monoxide. 
These abrupt transitions from low to high conversion 
steady-states had occurred early in the life of the 
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activity was stabilized, and thus remained very constant 
during the subsequent experimental work, the majority of 
which was performed at temperatures below 423 K. 

In order to determine the effects of both the severe 
reaction conditions and the initial pretreatment, it was 
decided to examine the behavior of a new twenty gram batch 
of the Pt-Pd converter catalyst under various operating 
conditions. Before being placed in the reactor, the batch 
OLecatalyst=was pretreated at 6257K™in flowing dry“air for 
approximately twelve hours. Thus, this catalyst was 
Subjected to the same pretreatment as the previous batches 
Ofecatalyst: 

After being placed in the reactor, the catalyst was 
Subjected to reaction conditions at 393 K which were similar 
to those under which the oscillations shown in Figure 7.13 
were obtained. It was found that the catalyst was very 
inactive and that only a single very low conversion steady- 
State existed. The reactor temperature was then raised to 
423 K, and it was found that the catalyst displayed oscilla- 
tory behavior which was qualitatively similar to that shown 
ineeioure 7.43 for Vohe previcus™bateh? Of catakyse aus com. 
The reaction was then driven to a low conversion steady- 
state by flooding the reactor with carbon monoxide, and 
following this the reactor was purged with oxygen. This 
type of treatment was very Similar to that which the first 


batch of catalyst had been subjected early in its life. 
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Once the oxygen purge was started it was observed that 
the carbon monoxide concentration initially decreased solely 
GvestoOuthe mixing effect; but that upon reaching a certain 
point it abruptly decreased virtually instantaneously to 
zero concentration. This abrupt transition was similar in 
many repects to that shown in Figure 7.14, and it was found 
that the catalyst activity had increased senate the 
abrupt transition. Upon setting the operating conditions 
back to those in effect before the reactor was flooded with 
CO, it was found that while oscillatory ‘behavior still 
occurred, it was qualitatively different in that the period, 
amplitude and shape of the peaks had changed. 

The reactor was then subjected to the above process 
Several times, and it was found that each abrupt transition 
had a smaller and smaller effect on the activity of the 
catalyst. Finally it was observed that the catalyst 
activity became constant and that it was no longer affected 
by further abrupt transitions from low to high conversion 
steady-states. Following this, the reactor temperature was 
lowered to 393 K, and it was found that when the feed condi- 
ELOnSewere set at those given in Figure /oisethau oseiiilae— 
tory behavior occurred which was qualitatively similar to 
The teeSHOWNs Ene 1LOUL Gama ics 

The above similarity of the oscillatory behavior was 
taken as an indicatiron that the activity of the second batch 
of Pt-Pd catalyst was now similar to that of the first 
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to stabilize the catalyst activity as long as the catalyst 


bulk temperature is kept below that at which the ‘reaction 


PEGE EeatMent eWas performed, 


TevmCoOnCluSsions 


From the experimental observations presented in the 


preceeding sections of this chapter, a number of comments 


can be made regarding the mechanistic details of oscillatory 


behavior. 


ie 


Two qualitatively different types of oscillations can 
occur during carbon monoxide oxidation. At low tempera- 
tures, relaxation oscillations which are related to the 
kinetic mechanism of the reaction predominate, while at 
high temperatures, small symmetrical oscillations caused 
by transport limitations occur. 

Chaotic behavior can occur when both types of oscilla- 
tions are occurring simultaneously due to the coupling 
of the two nonlinear oscillators. 

Since oscillations were observed at temperatures as low 
as 318 K, it appears that the reversible formation of 
surface oxides of platinum and/or the migration of 
oxygen into the bulk platinum must be ruled out as a 
possible cause of the relaxation oscillations. 

The presence of hydrocarbon impurities in the feed 
stream are eliminated as a cause of oscillatory behavior 


due to the existence of oscillations when ultra high 
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purity oxygen was used as the feed gas. Also, no water 
peaks were ever detected chromatographically, thus 
ruling out any significant accumulation and oxidation of 
hydrocarbons on the catalyst surface. 

Catalyst activity can be greatly increased by subjecting 
the catalyst to severe reaction conditions. It appears 
that this increase in the activity is a result of 
surface restructuring caused by high metal surface 
temperatures during periods of large instantaneous reac- 
tion rates. 

Since high instantaneous reaction rates occur whenever 
the carbon monoxide concentration oscillates from a high 
to a low level, then significant changes in temperature 
of the metal surface must occur whenever an oscillation 
occurs. However, if the catalyst has been ‘reaction 
pretreated! at a temperature which is higher than that 
at which the oscillations are occurring, then no changes 
invcatalyst activity will occurm due to the temperature 
Variactvons Guring the oscillations: 

It is believed that the oscillatory behavior observed at 
low temperatures is caused by an interplay between the 
rates of adsorption of carbon monoxide and oxygen, with 
thermal effects being produced by the large instanta- 
neous releases of energy during periods of abrupt tran- 
Sition from high"to Vow concentrations of "carbon 


monoxide. 
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8. A MODEL FOR OSCILLATORY BEHAVIOR FOR CO OXIDATION 


In the preceeding chapters reference has often been 
made to the many models which purport to explain oscillatory 
Dbehavecr, Geg. [Sn i7Saseo).9 Unfortunately; .thespredierions 
of these models have rarely been validated through detailed 
comparison with experimental data. In addition, many of 
these models contain simplifying assumptions, thus rendering 
their predictions suspect when viewed in the light of the 
development in Chapters 2 to 4. In this chapter, a model 
for the catalytic oxidation of carbon monoxide will be deve- 
loped which is free of those simplifying assumptions which 
cause violations of mass conservation, and some initial 
results from this model are compared to the experimental 


observations contained in the preceeding chapter. 


8.1 Model Development 

In the previous chapter it was stated that variations 
in the catalyst surface temperature were believed to play an 
important role in the observed oscillatory behavior. In 
Particular, the interaction of these temperacure variations 
with the rates of reaction and adsorption of oxygen and 
carbon monoxide was considered to be a key step in the reac- 
tion mechanism. Significant temperature rises during 


adsorption and reaction have been reported for supported 
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metal catalysts [8.7 - 8.9], and Dagonnier and co-workers 
[8.4, 8.5] have proposed a model that considers variations 
in surface temperature to explain oscillatory behavior. 
Their three ordinary differential equation model does not 
consider variations in bulk phase concentrations which, as 
has been shown in Chapter 7, are very significant, hence a 
more complete model for the catalytic oxidation of CO over a 
Supported metal catalyst was developed. This new model is 


based on the following assumptions: 


1. The experimental recycle reactor is assumed to behave as 
an ideal continuous stirred tank reactor. 

Zu (nem Leeagrrorthe reactor only contains ©O,, CO and inerts, 
1.e€. no CO, 1S present in the reactor feed. 

3. the réactor pressure’ is constant. 

4, The temperature of the bulk gas phase is constant and 
the feed stream enters at the reactor temperature. 

beeeehe Cemperaturesot the catalyst. support)\is identical co 
the bulk gas phase temperature. 

6. The temperature of the supported metal crystallites can 
vary and all heat transfer from the metal 1S approxi- 
mated by a single convection term. 

7, Thee physical*® properties of the metal@catalyst ane 
temperature invariant. 

8. The temperature of all metal crystallites is the same. 

9. The temperature of the metal catalyst is uniform 


throughout each crystallite. 
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10. Neither the adsorption of CO nor O, are activated. 

11. Oxygen dissociates upon adsorption on the surface. Each 
oxygen (0,) molecule thus requires two surface sites for 
adsorption, 

12. Carbon monoxide adsorbs in the linear form on the cata- 
lyst surface and thus each carbon monoxide molecule 
occupies a Single surface site. 

ise be conversion of carbon monoxide to carbonedioxide 
proceeds only through the surface reaction of carbon 
monoxide and oxygen surface species, i.e. solely 
Penomiun-rinshelwood kinetics: (8.10 — 16.12). 

14. The reaction is considered to be facile (structure 
insensitive). 

15. The surface area for adsorption is the same as the 
Surface area for heat transfer. 

VGemune thnerts and the CO, do not adsorb on the metal 
Surface. 

17. The surface sites are conserved. 

18. Neither the reaction activation energy nor the heats of 
adsorption or desorption are affected by the surface 
coverage. 

From the above assumptions, the kinetics of the reac- 
tion mechanism can be depicted as follows (S is a surface 


Site): 
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For this mechanism the use of the preceeding assump- 
tions leads to the development of a model which consists of 
Six ordinary differential equations, namely: three gas phase 
mass balances for CO, 0, and CO,; two surface phase species 
balances for COS and OS; and one overall energy balance for 


the metal catalyst. These six equations are written as 


follows: 
Balances! Onpgas phaser CO, OF and CO, 
VeolCOU—-O, CO} es OTCO}—=- mAb CONES]: + Ako pL COS] (8.4) 
ele 
V alo ] = On (G3% ie OLoen = ALO lode ae Res OS 12 (8.5) 
Ce 
VecgpCOsim=— = O| COs) teak [COST fos } (8.6) 
ake 
Balances On Surtace phase COS TandyO> 
BEeOtOS © Ss ZAK Oe Si 2 > -2Akos [OS Ie Ak ehCoOSs ll OSs (8¢7) 
oo 
Neo) COS. | ="Aky UCOGUS |) = Ak 4 GOS) se aki COSMGSs (8.8) 
(eas 5 


Energy balance on the metal catalyst 


mCp ane, = AAHE CeGEOl [Ss | = kee eos) mu Hrptks (Osis i2 
ene 
ei KOS ene Hrs ke ([COSThOS]} -eUAt a= T.)) (8.9) 
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These equations can be written in dimensionless form by 


defining the following dimensionless groups: 
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Inserting the dimensionless groups into Eqs. 8.4-8.9 gives: 
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8.2 Stability and Uniqueness 

An analytical investigation of the regions of stability 
and uniqueness for Eqs. 8.35-8.40 is obviously impossible 
due to the complexity of the system. Instead a numerical 
investigation was performed so that for any given set of 
parameters, the number and stability of the steady-states 


could be determined (see Appendix E for computer programs). 
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It was found, over large ranges of parameter values, that 
both steady-state multiplicity and instability were 
possible. Indeed, for some sets of parameter values as many 
as five steady-states were found to coexist, with varying 
numbers of the steady-states being stable or unstable 
depending on the parameter values. 

When large numbers of steady-states coexisted, the 
behavior of the system for a given set of initial conditions 
could only be determined by numerical integration (computer 
programs in Appendix E). From these integrations it has 
been found that this model can display complex dynamic and 
Steady-state behavior such as steady-state hysteresis, 
Simple periodic solutions, coexistence of stable limit 
cycles with stable steady-states, and multipeak oscillatory 
Henavrcone) "Arsoror incerest was the finding what; stor all 
the parameter values investigated, this model always 
possessed a stable low conversion steady-state (surface 


almost completely covered with CO). 


8.3 Numerical Integration Results 

In this section, selected results of numerical integra- 
f£ions Of Bas. 8.35-8.40 for several sets) Of parameter values 
are presented. The results are of a preliminary nature and 
should not be considered as an exhaustive delimiting of the 


behavior of the model. 
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Shown in Figures 8.1(a) and (b) are the results of 
numerical integrations of Eqs. 8.35-8.40 for a set of para- 
meter values for which five steady-states (four of which are 
unstable) existed. It is seen that the oscillations are 
qualitatively and quantitatively similar to the experimental 
BesuerSs=Showm lm Figure 7.3, as not sonly arerthe osctila- 
tions similar with respect to period and amplitude, but even 
the shapes of the oscillations agree in many respects. The 
major qualitative difference between the experimental and 
model oscillations occurs in the region of the carbon 
dioxide maximum. Experimentally the carbon dioxide concen- 
tration continued to slowly increase after the virtually 
instantaneous increase from the minimum, whereas the carbon 
dioxide concentration shown in Figure 8.1 decreases slowly 
after the instantaneous increase. However, the other 
regions of the experimental and model oscillations are in 
Substantial agreement. 

The Similarities between the model and experimental 
results are further emphasized in Figures 8.1(c) and (d) 
where the amplitude of the oscillations have been scaled so 
that the experimental order of magnitude difference in IR 
Sensitivity between CO and CO, 1s ‘duplrcated in theemodel 
results. 

If it is assumed that the mechanism postulated in 
deriving Eqs.7.8.35-9.40 15 an accurate representation of the 
physical reality, then the model can be used to determine 


the behavior of variables (temperature, concentration, etc.) 
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Figure 8.1 Oscillations Similar to Experimental Observations 
(parameters given in Table 8.1) 
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Table 8.1 Parameter Values for Figures 8.1 and 8.2 
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which were not measured experimentally. For example, the 
model predicts that the gas phase CO and CO, oscillations 
shown in Figure 8.1 are accompanied by catalyst temperature 
Occ Mreathonomot 520K itor U0 KhaGie0 K amplitude), and carbon 
monoxide surface coverage oscillations over the range 107° 
tomi0== {complete coverage = 1). It is in the determination 
of this type of information for variables which are diffi- 
cult (or impossible) to measure experimentally that a model 
can be particularly useful. 

The model predicts considerably different behavior than 
that shown in Figure 8.1 for other parameter values. Shown 
in Figure 8.2 are carbon monoxide oscillations for three 
different sets of parameter values. These oscillations 
illustrate the extremes which are possible, as either slow 
oscillations with a relatively broad maximum, or fast oscil- 
lations with sharp peaks are possible. A comparison of 
these oscillations with the experimental data shows that the 
Sscml lat Lonomin wrigunresed.,2\a) andmtb) aressimilanstouchose 
shown in Frgure 7.9. 

Besides the single peak oscillations shown in Figures 
Sf ana 8.2, this model can predict complex mulerpeak 
behavior as shown in Figure 8.3, where it 1S seen that two 
or three maxima per cycle behavior can occur. While the 
Goublet nature of the CO, oscillation is barely evident in 
Paogure 8.304), bt aS Very DLOnOUNced: InN PQULemO sb) Ramet 
BS) interesting eo note that the CO, oscillations shown in 


Figure 8.3(b) are qualitatively very similar to that shown 
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Figure 8.2 
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Extreme Forms of the CO Oscillations 
(parameters given in Table 8.1) 
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Figure 8.3 Multipeak Oscillations 


(parameters given in Table 8.2) 


Table 8.2 Parameter Values for Figure 8.3 
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in Figure 7.1(b). However, the complex three maxima per 
Gyaresoscillation shown in Figure 8.3(c) does not correspond 
to anything which was observed experimentally. These oscil- 
lations (three maxima per cycle) were produced by adjusting 
the parameters so that the model predicted large CO surface 
coverages, with correspondingly large amounts of CO desorbed 
from the surface during the high temperature parts of the 


oser lations, 


8.4 Conclusions 

The oscillations shown in Figures 8.1 to 8.3 are the 
results from a preliminary investigation of the model 
Composcdmotmnas. —o. 35-6.405. The anitial, results ci the 
model predictions are very promising since this model is 
able to reproduce some of the fine structure of the experi- 
mentally observed oscillations. However, considerably more 
model testing and experimental verification is required 
before definite conclusions about the validity of the model 
can be made. Nevertheless it can be definitely concluded 
that oscillatory behavior does not require multiple reaction 
mechanisms, heats of adsorption and activation energies that 
are functions of surface coverage, or the presence of unre- 


active adsorbed species. 
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8.5 Nomenclature 


Roman 


Surface area of metal 

heatmcapacitys Ofpthes catalytic metal 

Gas ‘phase concentration of CO 

gas phase concentration of CO, 

normalized surface phase concentration of adsorbed CO 
activation energies 

ratio of inlet to exit volumetric flow rates 

energy released upon adsorption or reaction 
pre-exponential miactors for =the: i-th reaction 

rate constants for the i-th forward reaction 

rate constants for the i-th reverse reaction 
dimensionless rate constants for the forward reactions 
dimensionless rate constants for the reverse reactions 
surface capacity (active sites per unit area) 

mass of catalytic metal 

gas phase concentration of oxygen 


normalized surface concentration of adsorbed oxygen 


‘pressure 


volumetric flow rate 

normalized concentration of active sites 
dimensionless time 

time 

dimensionless temperature 


metal crystallite temperature 
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tye temperature of the reactor feed 

U reactor heat transfer coefficient 

iV reactor volume 

Xx dimensionless concentration of carbon monoxide 
ry dimensionless concentration of oxygen 

Zz dimensionless concentration of carbon dioxide 
Greek 

®,,2 Gimensionless ratios of gas to surface capacitance 
on dimensionless heat transfer coefficient 

B; dimensionless heats of reaction and adsorption 
Yj dimensionless activation energies 

Subscript 


feed conditions 
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9, RECOMMENDATIONS 


Conclusions have been presented at the end of each 
chapter and will not be repeated here, but in the course of 
pursuing the research described in this thesis many 
additional areas of possible research have been encountered. 
While limitations of time, equipment and/or interest have 
resulted in these research areas not being fully pursued, it 
is realized that others may wish to do so. Listed below are 
recommendations for future experimental and theoretical work 
that will further the understanding of heterogeneous 


Gatabytrc reactions, 


1. Further investigations into the effects of simplifying 
assumptions are definitely needed. Of particular 
usefulness would be the development of criteria for 
determining a priori when specific simplifying assump- 
Eronstarée justified: 

2. General techniques should be developed for the 
determination of regions of steady-state uniqueness for 
systems of equations which obey conservation laws. 
Presently the regions of uniqueness can only be 
determined by various ad hoc procedures. 

3. The phenomenon of chaotic solutions for models of 
chemical systems needs further investigation. A 


systematic method of determining the regions (or even 
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the existence) of chaotic behavior in systems of 
equations would be a major achievement. At present, 
parameter values which give rise to chaotic behavior can 
only be wround by “chance”, P.e. through the -=numerical 
integration of equations for many sets of parameter 
values. 

During the numerical examination of oscillatory and 
chaotic behavior, it was found that the computer 
algorithms used to perform the integrations (Gear's and 
Runge-Kutta-Fehlberg) were inefficient. The development 
of specialized algorithms for integrating equations with 
oscillatory and chaotic solutions would be very useful. 
Additional reaction systems should be investigated to 
determine if oscillatory behavior over supported metal 
catalysts is restricted to oxidation reactions. 
Experimental measurements of surface concentration and 
temperature oscillations would be useful in further 
elucidating the reaction mechanism. For CO oxidation, 
in situ IR measurements could be used to determine the 
relationship between gas phase and surface phase 
concentration oscillations. 

The experimental data on the effects of operating 
conditions should be used to validate (or invalidate) 
the predictions of the many mogelsSewhicnspredice 
oscillatory behavior. 

Futher development and testing of the proposed six ODE 


model is warranted due to the initial promising results. 
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Appendix A Computer Programs for Chapters 3 and 4 


This appendix contains the source code of the program which 
was used to integrate the various systems of equations 
encountered in Chapters 3 and 4. This program was run on a 
Hewlett-Packard model 1000 21MX E-series computer, hence the 
code contains certain statements which are not part of 
Standard FORTRAN. These statements were necessary to 
perform the various file handling operations on the HP-1000. 
Also included in this appendix is the source code for the 
program which calculated the number and the stability of the 
steady-states for the three ordinary differential equation 
model derived in Chapter 4. Since this program uses 
suprouctnese:1rom the IMSL library, it had t@ be run on- the 
main university computer system (MTS - Amdahl), as the IMSL 


Subroutines are only available on that system. 
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PROGRAM IS USED TO INTEGRATE THE VARIOUS VERSIONS OF 
EIGENBERGER'S MODEL. SUBROUTINES F1, F2 AND F3 ARE USED 
TOPENTEGRATE, SYSTEMS 17 3 AND 2 (CHAP. 3) (RESPECTIVELY. 
SUBROUTINES F4, F5 AND F6 ARE USED FOR INTEGRATING 
SYSTEMS 1A, 3A AND. 2A REPECTIVELY, AND SUBROUTINE F7 IS 
USED IN THE INTEGRATION OF THE FULL MODEL (CHAPTER 4) 


i Qi @sG@@i@ @ Oe 


TN4 
PROGRAM EIGEN(,800), VERSIONS OF EIGENBERGER'S MODEL 
EXUBMPNAPEP ICP 2er ora F5 Po ry 
DIMENSION IDCB(144) ,NAME(3),YY(5),IBUF(100) ,NAMEI (3) 
*, SPARAM(6) , IDCBI (144) 
DOUBLE PRECISION Y(4),A,H,HMX,AER,RER,DA,YP(4) 
*, PARAM(6) ,DDA,HH,HHMX,ASV,YSV(4) ,DABS,DSORT 
COMMON PARAM 
GAVE ADTACH 
LUNW=LOGLU(ISIS) 
LUNR=LUNW 
CALL IFMTE(LUNW) 
WRITE(LUNW, 800) 
800 FORMAT(" ENTER NAME OF INPUT FILE, AND LU FOR OUTPUT") 
READ(LUNR,290) NAMEI,LUNO 
290 FORMAGUGAZ. bho) 
CALL OPEN(IDCBI,IERR,NAMEI,0,1) 
CALL READF(IDCBI,IERR,1BUF) 
GAT GaCODE 
READ(IBUF,*) A,DA,AER,RER 
CALL READF(IDCBI,IERR,IBUF) 
GAGE EGODE 
READ(IBUF,*) N,NF,NP,NINT,IPLOT 
NOV=5 
NIV=4 
CALL READF(IDCBI,IERR,IBUF) 
GAG ign CODE 
BEAD CEBU Ss mci Gb yl i eNiV) 
CALL READF(IDCBI,IERR,IBUF) 
GALIn CODE 
READ(IBUF,*) (PARAM(I),1=1,NP) 
WRITE(LUNO, 901) 
CALL DFILE(IDCBI,IDCB,NAME) 
'GALEECORE 
WRITE(IBUF,350) NINT 
350 FORMAT (314) 
CALL WRITF(IDCB,IERR,IBUF,2) 
GALImMCODE 
WRITE(IBUF,350) N,NF,NP 
CALL WRITF(IDCB,IERR,IBUF,6) 
DOs (Uist = | phe 
10 SPARAM(I )=SNGL(PARAM(1) ) 


GALGECODE 
WRITE(IBUF,361) (SPARAM(I),I1=1,NP) 
361 FORMAT (6E12.5) 


CALL WRITF(IDCB,IERR,IBUF,6*NP) 


‘ a : 
(<, WaT aNe? 
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DDA=DA 

HH=DDA/4.0D0 

HHMX=DDA/2.0D0 

YY C1) =SNGL CA} 

DOMT29i=1 NIV 

Vy Ul tiie SNGE (Yn) 

CHEUACODE 

WOU E BURY GG Tim YY.) ) Lec NOV) 
CALL WRITF(IDCB,IERR,IBUF, 6*NOV) 
NPINT=1 

ICOUNT=0 

DOAS ORY =| NINT 

I DUM= 1 

CONTINUE 

DO 26 KK=1,IDUM 


l=1 

SV ety Y CL) 
O20 30 7405756 76047 10) NE 

CALL RKF(A,Y,N,F1,DDA,HH,HHMX,AER, RER, IFLAG) 

¥ (4) =0.0D0 

¥(3)=1.D0-Y(1)-Y(2) 

COTO. 700 

CALL RKF(A,Y,N,F2,DDA,HH,HHMX,AER,RER, IFLAG) 

Velo 70D0 

Wee = ep ay Cl) —V 63) 

GORTO R790 

CALL RKF(A,Y,N,F3,DDA,HH,HHMX,AER,RER,IFLAG) 

Y¥(4)=0.0D0 

Via vad pO-¥ (2)-Y(3) 

CGOITSTTG0 

CALL RKF(A,Y,N,F4,DDA,HH,HHMX,AER,RER,IFLAG) 

¥(3)=1.D0-Y(1)-Y(2) 

TR(PARAM(5)7LE. 120D-10)°GO TO 745 

Y(3)=(-1.D0+DSORT(1.D0+8.D0*PARAM(5)%*(1.D0-Y(1) 

#=¥(2))))7(4.D0*#PARAM(5) ) 

Y(4)=PARAM(5)*Y(3)*Y(3) 

GOGTO 3790 

CALL, REE CA, Y> 

weer =). DOK v4 

Y(4)=PARAM(5) 

‘GO LTO 17.90 

CALL RKF(A,Y,N,F6,DDA,HH,HHMX,AER, RER, IFLAG) 

V(t j= Tt DO=Y 02) Sy (3) -2.DOX*XPARAM(5) 4¥ (3) 4¥ (3) 

Y(4)=PARAM(5)*Y(3)#Y(3) 

GO UTO 1790 

CALL RKF(A,Y,N,F7,DDA,HH,HHMX,AER,RER,IFLAG) 

¥(3)=1.D0-Y(1)-Y(2)-2.D0*v(4) 

LE (LF EAGUNE. 18 GO eTO 1010 

YY(1)=SNGL(A) 

bOe25 Vda, NL. 

Vvidd+ 1) =SNGL(Y (Jd) ) 

IF(KK.EQ.IDUM) GO TO 26 

CALURCODE 


F5,DDA,HH,KHMX,AER,RER, IFLAG) 
Y¥(3)-2.D0*PARAM(5)*Y(3)*Y(3) 
( 
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26 


16 


27 


900 
SO 
9US 
30 


G10 


100 
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Wie BUR oq meee) a= | NOV) 
CALL WRITF(IDCB,IERR,IBUF,6%*NOV) 
We EE OCUNO | O0'S CV Yat) r= 1 NOV.) 
NPINT=NPINT+1 

CONTINUE 

BeCLUUMeNE 21.) GO ro med. 
PLAS GUY OV Mov Woy iGloninur sO. 10D0)) GO TO 27 
A=ASV 

DO la t= NIV: 

Veale vsv C1.) 

I DUM=10 

BOA=DA AIO. ODO 

HH=DDA/4.0D0 

HHMX=DDA/2.0D0 

GORTO™? 2 

CONTINUE 

I COUNT=ICOUNT+1 

DDA=DA 

HH=DDA/4.0D0 

HHMX=DDA/2.0D0 

Be @lPCOUNT 4IPLOT+IPLOT. NEx PCOUNTOY GOUTO 30 

CALMECODE 

WRITE(IBUF, 361) (YY¥(1),1=1,NOV) 

CALL WRITF(IDCB,IERR,1IBUF,6*NOV) 

Were MUNG, 90>) VY CIl 3 G= NOV) 

NPINT=NPINT+1 

BORMAGO> Xu MED reek Aen Ck Bhat ake Ke) 

Megs 5 en UME ee le A A UU Kon Xx ae XBR) 


RORMAT (6 0Ue E4000, 77 k2 O4P15. 6, 1k) 
CONTINUE 

WRITE(LUNO,910) NPINT 

BORMAT Cm ee Or POTNTS =" 1 5) 


GALE EChOSE (1 DCB) 

CALL OPEN(IDCB,1IERR,NAME,2,1) 
CALL RWNDF(IDCB) 

CREE CODE 

WRITE(IBUF,350) NPINT 

Gro eAR TE (TOCE ATERR sPBUPF 32) 
GChbE CLOSE t(rDCB) 


STOP 

CONTINUE 

WPT TE (LUNO, 915) (§EUAG AyiY tll) llth 
BORMAP UD pum iG Glis, 6a) 


Chie CLOSECEDGS ) 

CALL OPEN(IDCB,IERR,NAME,2,1) 
CALL RWNDF(IDCB) 

GARD GOD 

WRITE(IBUF,350) NPINT 

CAL WeLEE IDC o eV LERE, LBUR wey 
CALEECEOS= GELDER) 

STOP 

END 

SUBROUTINE F1(T,Y,YP),DAX/DT, DBX/DT (ALL ASSUMPTIONS) 
DOUBLE PRECTSION Y(3) /YPCU3Q) 7 
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DOUBLE PRECISION AX,BX,DAX,DBX,D,E,F,H,G,GR 

COMMON D,E,F,H,G,GR 

AR= Vil) 

BX=Y (2) 

DAX=D* (1.D0-AX-BX)-E*xAX-F*AX*AX*(1,D0-AX-BX) **2 
DBX=H*(1.D0-BX-AX)-BX 

VERA Saw 

YP(2)=DBX 

Vesa 0 DU 

RETURN 

END 

SUBROUTINE F2(T!]V)VP) ) DAX/DT, DX/DT: (ALL ASSUMPTIONS) 
BOUBLE PRECISION-Y(3) -YP(3)57T 

DOUBLE PRECISION AX,X,BX,DX,DBX,D,E,F,H,G,GR 

COMMON D,E,F,H,G,GR 

Ax=V (1) 

Xeay (3) 

DAX=D*X-E*AX-F*AX#AX*X#X 
DX=-D*X+E*AX-H*¥X+1.D0-X-AX+2.DO#F*AX*AX¥X#X 
YP(1)=DAX 

Ves) =r DO 

YPC3)=Dx 

RETURN 

END 

SUBROUTINE F3(T,Y,YP),DBX/DT, DX/DT (ALL ASSUMPTIONS) 
POURLEEPRECTSLON ey (3) YPC3) 3 

DOUBLE PRECISION X,BX,DX,DBX,D,E,F,H,G,GR 

COMMON D,E,F,H,G,GR 

Vales ®) 

Bray (2) 

DBX=H*X-BX 
DX=-D*xX+Ex*(1.D0-X-BX)-H*X+BX+2,D0*«F*(1.D0-X-BX) **2*X*X 
YP(3)=DX 

ve (2) =DBA 

VEG 0%, 0D0 

RETURN 

END 

SUBROUTINE F4(T,Y,YP),DAX/DT, DBX/DT EQUIL. ASSUMPTION 
DOUBLE PRECISION Y¥(C3)-YPt3) 37 

DOUBLE PRECISION AX,%,BX DAX DX, DBS, DE FH GeyGR. DSORT 
COMMON D,E,F,H,G,GR 

AX=VC 1) 

BRaViee) 

X=1.D0-AX-BxX 

Rea (se hee) OD i GO. Ort 0 
X=(-1.D0+DSORT(1.D0+8.D0*G*(1.D0-AX-BX) ))/(4.D0%G) 
DAX=D*¥X-E*xAX-FXAX*XAX*K*X 

DBX=H*X-BX 

YP(1)=DAX 

YE CZ)=DBx 

YEOoJ=0.0D0 

RETURN 

END 

SUBROUTINE F5(T,Y,YP),DAX/DT, DX/DT EQUIL. ASSUMPTION 
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DOUBLE’ PRECISION Y(3) ,YP(3),T 

DOUBLE PRECISION AX,BX,X,DAX,DBX,DX,D,E,F,H,G,GR 
COMMON D,E,F,H,G,GR 

Asay (M1) 

a Sy) 

B= 1-D0-AR—X-2 DOFGHEREX 
DAX=D*X-E*AX-F*AX*AX¥X#X 

DX=-D&*X+E* Y-H*X+BX+2 .DOXFXAAKXAK#X#X 

YOU?) =DAs 

¥PU2) =0 s0D0 

YeGs)=Dx 

RETURN 

END 

SUBROUTINE F6(T,Y,YP),DBX/DT, DX/DT EQUIL. ASSUMPTION 
DOUBLE, PRECISION Y(3) -YPC3)>T 

DOUBLE PRECISION AX,X,BX,DX,DBX,D,E,F,H,G,GR 
COMMON D,E,F,H,G,GR 

Rays) 

BX=Y(2) 

AX=1-BX-X-2.D0*G*X*xX 

DBX=H*X-BX 
DX=-D*X+E*AX-H*¥X+BX+2.DO*FXAX*XAXXX*X 
YP( 3) =Dx 

YP(2)=DBX 

YP(1)=0.0D0 

RETURN 

END 

SURROUDINE aha lay “Vek DES DEBk7 DT wAND DBXX/DT 
DOUBLE PRECISION Y¥(4),YP(4),T 

DOUBLE PRECISION AX,BX,XBX,X,DAX,DBX,DXBX 
DOUBLE PRECISION D,E,F,H,G,GR 

COMMON D,E,F,H,G,GR 

ARH) 

BX=Y(2) 

XBX=Y(4) 

X=1.DO0-AX-BX-2.D0*XBX 
DAX=D*X-E*xAX-F*AX*AK*XBX/G 

DBX=H*X-BX 
DXBX=GxGREX*¥X-GR*¥XBX-FXAXXAX*¥XBX/(2.D0%G) 
YP(1)=DAX 

WZ) =DESx 

YP(4)=DXBX 

me oy —0 2000 

RETURN 

END 

SUBROUTINE DFILE(IDCBI,IDCB,NAME) ,OUTPUT FILE CREATION 
DIMENSION IDCB(144) ,NAME(3),ISIZE(2),IDCBI(144) 
CALL READF(IDCBI,IERR,NAME) 

CALL CLOSHUPDEBT) 

Ol ZE VL) =45 

PIYen=s 

PSECU 

TCR=135 

CALL OPEN(IDCB,IERR,NAME,0,1) 


CLOL@raw@1@) 


Oe 


—J 


CALL CLOSE(IDCB) 

CALL PURGE(IDCB,IERR,NAME, 1) 

CALL CREAT(IDCB,IERR,NAME,ISIZE,ITYPE,ISECU,ICR) 
CALL OPEN(IDCB,IERR,NAME,0,ISECU) 

RETURN 

END 

SUBROUTINE RKF(A,Y,N,F,DA,H,HMX,ABSER,RELER,IFLAG, 
*,NPINT,IDCB), VARIABLE STEP INTEGRATION OF ODE'S 


THIS CODE INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY 
DIFFERENTIAL EQUATIONS BY RUNGE-KUTTA-FEHLBERG METHOD 
WITH AUTOMATIC ESTIMATION OF LOCAL ERROR AND STEP SIZE 
ADJUSTMENT - SHAMPINE AND ALLEN (REF. 5.16) 


DOUBLE PRECISION Y(N),YTEMP(20) ,TEMP(20),R(20) 
DOUBEEMPRECTSLONSR PC2 0) 7K 2 G20))4eK 3020) PReE(20) 7 K5(20) 
*,K6(20),U,RELER,ABSER,HMX,A,B,DA,HMAX,H,HKEEP,ARG 
*,RATIO,DABS,DMAX1,DMIN1,DSIGN,DSORT,TR 

REAL DOUT(20) 

INTEGER IBUF(100) ,IDCB(144) 

US. 90D t1 

Pr -CRELER. LPO 20D0.O0R.ABSER.LT.0.0D0) GO TO 18 

IF (RELER+ABSER.EQ.0.0D0) GO TO 18 
PEXHMX GE J0S0D0), "Go. TO 18 

B=A+DA 

PE CDABS (DA) LETS .DO*U#DMAX1 (DABS(A) ,DABS(B)))GO TO 19 
HMAX=DMIN1(HMX,DABS(DA) ) 
IF(DABS(H).LE.13.0D0*U*DABS(A)) H=HMAX 

IADJS=0 

H=DSIGN(DMIN1(DABS(H) ,HMAX) ,DA) 

IF (DABS(B-A) .GT.1.25D0#DABS(H)) GO TO 4 

HKEEP=H 

PADGS=1 

H=B-A 

Geiger GAs Y AKT) 

CONTINUE 

DGA6m = 1 N 

YTEMP(1I)=Y(1)+0.25D0*H#K1(1) 

ARG=A+0.25D0*H 

CALL F(ARG,YTEMP,K2) 

DOS BEN 
YPEMP (yal ptHe( Ki1G1 ) #3 2D0 432: DOS KSI US D07 32. D0) 
ARG=A+H*(3.D0/8.D0) 

GCARLVIE.CARG-YTEMP K3) 

DOy Gre =a, N 
YTEMP(1I)=Y(1)+H*(K1(1)*(1932.D0/2197.D0) 

= Kidie ae 720020072197. D0) +KS (Lj 29 GR aUy 27:97 3G) e) 
ARG=A+tH*(12.D0/13.D0) 

CALL F(ARG,YTEMP,K4) 

Oma =i 
YTEMP(1I)=Y(1)+H*(K1(1)*(439.D0/216.D0)-8.D0*K2(1) 
*+K3(1)*(3680.D0/513.D0)-K4(1)*(845.D0/4104.D0) ) 
ARG=A+H 

CALL F(ARG,YTEMP,K5) 


1 


a5 
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Omir t=d oN 

VIEMP( 1) =¥ OL) thet-K 1 (1) ¥(8)D0/277D0) +2 sD0*K2(1) 
*-K3(1)*(3544.D0/2565.D0)+K4(1)#(1859.D0/4104.D0) 
Poh oes aie D040. DON) 

ARG=A+0.5D0*H 

CALL F(ARG, YTEMP,K6) 

DOM its = 0 RN 
TEMP(I)=K1(1)*(25.D0/216.D0)+K3(1)*(1408.D0/2565.D0) 
*+K4(1)*(2197.D0/4104.D0)-0. 2D0*#K5(1) 
YTEMP(1I)=Y(1)+H*TEMP(I) 

DOG h=4-N 

RGI=K (1) / 360eD0-K3 (1) *(128.D0/4275.D0)-K4-(7 } 
ee 1972D0/75240.D0)+K5(1 ) 750). DO+KG6(1 )*(2.D0/55.D0) 
RATIO=0.0D0 

DORIS ras TEN 

TR=DABS(R(1))/(RELER*DABS (YTEMP(I))+ABSER) 
RATIO=DMAX1(RATIO,TR) 

PRIGRATIO.GT. 17D0) GO TO 15 

DO 40 I=1,N 

DEAT) FR SCID RG .0: ODOR GOTO “40 
Pho BO 202 0NG) EGO TORet 

DOUT(1)=SNGL (A) 

DG seo I= N 

DOUT(Gd+1)=SNGE (CY (Jd) ) 

M=N+1 

TF (N. EOS 2) M=N+2 

CALE CODE 

WRITE(IBUF,900) (DOUT(JJ) ,JJ=1,M) 

FORMAT(11E£12.5) 

CALL WRITF(IDCB,IERR,IBUF, 6M) 

NPINT=NPINT+1 

GO TO 45 

CONTINUE 

CONTINUE 

DOw141=i19N 

Y(1I)=YTEMP(1) 

A=A+H 

LP CLADISSEOG SIN GOSTORIE 
RATIO=DMAX1(RATIO,6.5536D-04) 
RATIO=DMIN1(RATIO,4096.D0) 
H=0.8D0*H/DSORT(DSORT(RATIO) ) 
IF (DABS (H).LE.13.D0#U*DABS (A) 
LECRATUOSRS D0 RGOLTSG IS 
PAD S=0 

GORTG?5 

IFLAG=1 

H=HKEEP 

RETURN 

IFLAG=3 

RETURN 

LEDAG=4 

RETURN 

END 


peereuar cee! 
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THIS PROGRAM IS USED TO DETERMINE THE NUMBER AND 
STABIEITY OF THE STEADY-STATES FOR THE GENERAL 
THREE ODE MODEL DERIVED FROM EIGENBERGERS MECHANISM 


EXTERNAL F 

DOUBLE PRECISION Y(1),WA(3),PAR(6),EPS,DSORT,F 
DOUBLE PRECISION AX,BX,XBX,X,DAX,DBX,DXBX 
DOUBLE PRECISION FAC1,FAC2,FAC3 

DOUBLE PRECISION RJ(3,3),WK(3) 

COMPLEX#16 W(3),2(3,3) 


READ(5,790) NRUN 
790 FORMAT(I2) 
DO 100 KK=1,NRUN 


#35 
794 


136 


~ 


800 
900 


* 


810 
905 


o 1 


READ(5,795) ITMAX,NSIG,EPS 
FORMAT(213,D10.5) 
WRITE(6,794) 
RORMAT (| UR) (77 ) 
WRITE(6,796) ITMAX,NSIG,EPS 


BeeMAT GT | Oe MA Nera TRRATRONS a= 6 i] Sees eee STG. 7 
eG Cee deen ee MA XP UNGTS aVALUP ==> DiseG, (7) 
READ (5, 200) (PARC), b=1,6) 

FORMAT(6D10.5) 
WRITE(6,900) PAR 

PORMAT (” OSE TO! ek tae Se) LBC -dcnag a 

at Ole lam 52k See Ds Gree Oe TALE PAR =m Dis) 647 
TO ARG cme Dao Gere 

GE ey Gee =P Ah i/o TO Kh 5% Sys Waid ge) 
READ(5,810) Y(1) 

FORMAT (2D10.5) 
WRITE €6,905) ¥Y 

FORMAT oar, VIN TC LALALVARUB@Oles Go Xsr =n aD fol tos or) 
N= 1 

CALL ZSYSTM(F,EPS,NSIG,N,Y,ITMAX,WA,PAR,IER) 
WRITE(6,910) ITMAX,IER 

BORMAT 04) oO Re PCE RATIONS = Cel oy okie. © Re eee 


REX=Y( 1) 

FAC1=2.D0*PAR(1)*PAR(5)/(1.DO0+PAR(2) ) 
FAC2=(PAR(1)**2)*PAR(5)/(2.0D0*(PAR(4) **3) ) 

FACS=8. 0D0eFAC Ie{ 1.D0-2% DUSXEXG XB 
AX=(-FAC2-FAC1+DSORT( (FAC1+FAC2) **2+FAC3))/(4.D0*XBX) 
BX=(1,.D0-AX-2.D0*XBX) *PAR(2)/(1.D0+PAR(2) ) 
X=(1.DO0-AX-BX-2.D0#XBX) 
DAX=PAR(3)*PAR(1)*#X-PAR(3)*PAR(1)**2*AX/(4.D0x 


*PAR(4)**3)-PAR(3) *AX*¥AX*XBX/PAR(5) 


DBX=PAR(2)*X-BX 
DXBX=PAR(5)*PAR(6) *X*X-PAR(6) *XBX-PAR(3) ¥AX*AX*XBX/ 


*(2.D0*PAR(5) ) 


95:0 
* 


960 
*k 


WRITE(6,950) AX,BX,XBX,X 


FORMAT. UPTO PAxe  =)). Dimon sain? One enX mer ies rccs 
va Oe eX Sa =o ay yal ee Ey MR Serge 
WRITE(6,960) DAX,DBX,DXBX 

PORMAT UD ROpeeeDAS = DAS Gye Os OB Xue Dil sub, 
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RJ(1,1)=-PAR(3)*(PAR(1)+PAR(1)**2/(4.DO*PAR(4) **3) 
*+2,DO*AX*XBX/PAR(5)) 
RJ(1,2)=-PAR(3)*PAR(1) 


RJ(1,3)=-PAR(3)*(2.D0*PAR(1)+AX*AX/PAR(5) ) 

Red 02% 1) ==PAR(2) 

Rd(2)2)=-PAR( 21-100 

RJ(2,3)=-2.D0*PAR(2) 

RJ(3,1)=-2.D0*PAR(5) *PAR(6)*X-PAR(3) *AX*XBX/PAR(5) 
RJ(3,2)=-2.D0*PAR(5) *PAR(6) #X 

RJ(3,3)=-4.D0*PAR(5) *PAR(6)*X-PAR(6)-PAR(3) ¥AX*AX/ 


Ca DO*PAR(S)) 
GABU@ETGRE CRI, G5 oy0e Ww, 27 oe WK, LER? 
WRETEC6ONITO)) Go Wi Jaomodade 3) 
970 BORMAT.CTiT 1 ORT GENVALUES Hie dlls = Ds .6) + ly Dig. om led 
100 CONTINUE 
STOP 
END 
FUNCTION F(Y,K, PAR) 
DOUBLE PRECISION F,Y(1),PAR(1),DSQORT 
DOUBLE PREGCLSION X7AX, BRO XBR PACT FACZ FACS 
Pee 7) PLES, Ob0gm VG k= 120D—42 
ie ey ate 
FAC1=2.0D0*PAR(1)*PAR(5)/(1.D0Q+PAR(2) ) 
FAC2=PAR(1)**2*PAR(5)/(2.0D0*PAR(4) #*3) 
FAC3=8.DO*#FAC1*(1.D0-2.D0*#XBX) *XBX 
AX=(-PAC2-FAC1+DSORT( (FAC2+FAC1)**2+FAC3))/(4.D0*XBX) 
ht (1.,DO-AX-2. PD OREEM) RAR 2 AE) DOSPARA2 
X=1.D0-AX-BX-2.D0*«XBX 
Fe See tec X*X-PAR(6)*XBX-PAR(3) xAX*AX*XBX/ 
*(2.DO0*xPAR(5)) 
RETURN 
END 


Ac 


Appendix B Computer Programs for Chapter 5 


This appendix contains the source code of the computer 
programs which were used in Chapter 5. The first program, 
DIM3, was used to integrate the system of three ordinary 
differential equations using the RKF and Gear's methods. 
The second program, AUTOC, was used to calculate the 
autocorrelation coefficients. Both of these programs were 
run on the HP-1000 minicomputer, hence the deviations from 
Standard FORTRAN. Also included in this appendix is the 
source code for the program which was used to determine the 
regions of uniqueness and instability for the three ordinary 
differential equation model presented in Chapter 5. This 
program was run on the main university computer system (MTS 


- Amdahl) due to its usage of the IMSL Subroutines. 
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eC 
C THIS PROGRAM IS USED TO INTEGRATE THE SYSTEM OF THREE 
C ORDINARY DIFFERENTIAL EQUATIONS PRESENTED IN CHAPTER 5. 
C THE INTEGRATION CAN BE PERFORMED USING EITHER GEAR'S 
C METHOD OR THE RKF METHOD. 
© 
FTN4,X 
BLOCK DATA NAME, NAMED BLOCK DATA AREAS 
DOUBLE PRECISION DA1,DA2,B1,B2,GAM1,GAM2,ALPHA,DA 
COMMON /PARAM/ DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 
COMMON /CT/NPINT,IPINT,NDLAY,IMAX,IBUF(100) ,IDCBO( 144) 
END 
PROGRAM DIM3(,31005), THREE ODE MODEL (CHAOS) 
EXTERNAL F2,F3,PDERV 
DIMENSION NAMEO(3),X(4),NAMEI(3),IDCBI (144) 
*, IBUFA(50, I RBUF (33) 
DOUBLE PRECISION Y(3),A,H,HMX,AER,RER,DA 
DOUBLE PRECISION DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 
DOUBLE PRECISION TPT,YCO,ZHC,G(8,3) 
DOUBLE PRECISION GSAVE(39) ,GMAX(3),GERR(3),GPW(18) 
COMMON /PARAM/ DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 
COMMON /CT/NPINT,IPINT,NDLAY,IMAX,IBUF(100),IDCBO(144) 
EQUIVALENCE (NAMEI(1),IRBUF(2)),(ICRI,IRBUF(6)) 
EQUIVALENCE (ICRO,IRBUF(10)),(LUNO,IRBUF(14)) 
CALL DTACH 
LUNW=LOGLU(ISIS) 
LUNR=LUNW 
CALL IFMTE(LUNW) 
GAD AGERST (LBUPA ,= 17 ;1L0G) 
CALL PARSE(IBUFA,ILOG,IRBUF) 
EEGEUOG GT. 0) 4GO) TO, 600 
WRITE(LUNW, 800) 
800 FORMAT(" NAME OF INPUT DATA FILE ") 
READ(LUNR,290) NAMEI 
290 FORMAT (3A2,1X%,12) 
WRITE(LUNW,801) 
801 FORMAT("ENTER INPUT CRN, OUTPUT CRN, AND OUTPUT LUN") 
READ(LUNR,*) ICRI,ICRO,LUNO 
600 CALL OPEN(IDCBI,IERR,NAMEI,0,1,ICRI) 


CALL READF(IDCBI,1IERR,I1BUF) 
GALE, GODE 

- READ(IBUF,%*) N,INTG 
CALL READF(IDCBI ;,IERR,IBUF) 
CALL CODE 
READ(IBUF,*) A,DA,H,HMX,AER,RER 
CALL READF(IDCBI,IERR,IBUF) 
GALL ECODE 
PEAD( GBUR MH) “DPT .YCO, ZHG 
CALL READF(IDCBI,IERR,IBUF) 
GAELIGODE 
READ(IBUF,*) DA1,B1,GAM1,ALPHA 
EE NBEO 2 MGOsTOsT 
CALL READF(IDCBI,IERR,IBUF) 
CARI. GODE 


J 
iM 


er if voirvers ae  seeainclavent ‘ot eeu: 
ops wy CPWRSAG SOT Tadd CAL Rie 
2° tag) WRT YS aati Commute. Ut HAD Bee 
s4) ran de, 


Seo ARAC SOC COMA PRM AR Gas 
ome Kh KAD te. 10, At ee MO P23 e 
aca Mam Meo, SA, &. CAC a SAR eee 2 
CAGE) Vg Cae. Ode, Say Ce es) as 


oO) nt dee ee, igo Vii 
vee gt Ce Fara 
meet OPN )e MAR Wt 1ReNIS | * 
ly “att. TehA Se Ly wi : 
“ophe Sa A a BEDE BOTS Te eee 
i 8 OG: a "i 
Fi \! QP war ol23 DR SCS 7 7 
a b ap 14 OY Ay A LOSS ey ae <3 - 
"oe s, i. Bho, A AMAR Re 
a 7t ce pi ( WE (GUN TD MO - 
Qe. cr) rave) SRS RAvEDOS 7 
way « y § HUMML OF DT) BOIAVICRs rr 
PATS Jd7A2 -_ 
2.6 112 SGlenr 
Wie {aid oo 
a 
7 


by 


pers RIM Late 
KK, Sts AUT YDS aaa 


is ee ee 

P COS Wala at tae a 

fica PAG UAT: OO eA oe 
aran. 290% ABTA 7 

di. CAR) TAMAS ges 7a 

| Fis rier aT ie ' 7 
Wats | Per ep). reser. aes Rete" pao: 5) ee 
Wt ONS, Sie |e a ot 
| Cine, ee gs eS oe 

aren md ta Sah 3 


is 
_» 
mn BY ae 


a, . : io 


baat, tioat 4 +H 


READ(IBUF,%*) DA2,B2,GAM2 
1 CONTINUE 
CALL READF(IDCBI,1IERR,I1BUF) 
CABLE eCODE 
READ(IBUF,*) NINT,IPLOT,NDLAY 
WRITE(LUNO,330) A,DA,H,HMX,AER,RER 
WRITE(LUNO,330) TPT,YCO,DA1,B1,GAM1,ALPHA 
DF (NE BOs) GWRITE(LUNO #330) ZHC)DA2,B2,GAM2 
230 FORMAT(12(1PD13.5)) 
WRITE(LUNO,345) NINT,IPLOT,N,INTG,NDLAY 
345 FORMAT (6(5X,14)) 
a Es 
mez =¥ EC 
VCS) =ZHG 
X(1)=SNGL(A) 
Re = SNGE CY (1) 
RAS =SNGU(Y (2) 
KA = SNGLUY to) 
CALL DFILE(IDCBI,ICRI,IDCBO,NAMEO,ICRO) 
GALE CODE 
WRITE(IBUF,350)NINT 
350 FORMAT (16) 
CALL WRITF(IDCBO,IERR,IBUF,3) 
SDA1=SNGL(DA1) 
SDA2=SNGL(DA2) 
SB1=SNGL(B1) 
SB2=SNGL(B2) 
SGAM1=SNGL(GAM1) 
SGAM2=SNGL(GAM2 ) 
SALPHA=SNGL (ALPHA) 
CALL CODE 
WRITE(IBUF,360) SDA1,SB1,SGAM1,SALPHA 
360 FORMAT (12(1PE13.5)) 
CALE -WRITF(IDCBO,IERR,DBUF, 26) 
[rin EO 2 GO-TO?! 2 
CALL GODE 
WRITE(IBUF,360) SDA2,SB2,SGAM2 
CALL WRITF(IDCBO,IERR,IBUF, 20) 
2 CALL CODE 
WRITE(IBUF,360) A,DA,H,HMX,AER,RER 
CALL WRITF(IDCBO,IERR,IBUF,40) 
IF(NDLAY.GT.0) GO TO 5 
CALE ECODE 
WRITE CIBUR,S60) -(X(1),1=1,4) 
CADL WRITE (IDCBO,1IERR, BUF; 26) 


5 NPINT=1 
IF(NDLAY.GT.0) NPINT=0 
DOMAOL=A EN 
GMAX (1 )=0201D0 
7 Gar ee YT) 
JS=0 


DO 10 IPINT=1,NINT 
PFUNSEO, 2.ANDeINTG.EODO0) 
*CALL RKF(A,Y,N,F2,DA,H,HMX,AER, RER, IFLAG) 


Rat 


fi GRE Onn 133% 


su 


on cy (va 


e472 yo uma 
. @' wre 
is, a , 
a] nel 


% 3 Fick rst YP Ase 5 
, i ; =) A 
; aes ¢ “i 

7 ° ~~ 


LO 0 


10 


Ste) 


100 
Sts 


TP (NEO. 30AND.INTG.EO.0) 
*CALL RKF(A,Y,N,F3,DA,H,HMX,AER, RER, I FLAG) 
IF(INTG.EQ.1)CALL SGEAR(N,A,DA,G,GSAVE,H,HMX 
*,RER,1,GMAX,GERR,IFLAG,JS,6,GPW,PDERV) 
IF(INTG.EQ.0) GO TO 9 
DOCS ois In 
Vet = G Cais) 
TEQIFLAG SNE. 1) oGO TOS100 
PRURRINT.GT NDLAY) GO £0 10 
xi) =SNGL(A) 
MG) = SNGE CY (1) 
) 
) 


ROS) =SNGE Gy (2 

X(4)=SNGL(Y(3 
IF(IPINT/IPLOT*IPLOT.NE.IPINT) GO TO 10 
TR CLMAX (EQ 1) (GO TO +10 

CALL CODE 

WRETECDRUR 360)" (X(T) el =7)4) 
GALBDIWRETE(1DCBO; TERR, 1 BUF, 26) 


) 
) 
) 


Mee TDaVvEUNOvsS Oy mA YC) ) YY C2) 4 (3) 
NPINT=NPINT+1 

CONTINUE 

WRITE(LUNO,910) NPINT 

FORMAT ( Pao sear Gpo LNs) ES) 


CAL ELCLOSECIDCBO) 

CALL OPEN(IDCBO,IERR,NAMEO,2,1,ICRO) 
CALL RWNDF(IDCBO) 

GCALBECODE 

WRITE(IBUF,350) NPINT 

CALL WRITF(IDCBO,IERR,IBUF,3) 

CALL CLOSE CIDCBO) 


STOP 

CONTINUE 

WRETE (LUNG 190 5S) “1 BDAG PAY Gt) Va) ey (3) 
PORMAT GUN 13 /4(D15.6, 14)) 


CALM FEbOSE(GTDCBO) 

CALL OPEN(IDCBO,IERR,NAMEO,2,1,ICRO) 
CALL RWNDF(IDCBO) 

CAUET CODE 

WRITE(IBUF,350) NPINT 

CABDUEWRITE (PDGBO IERRVI BUF. 3) 

GALE CLOSE (1IDGBO) 

STOP 


' END 


SUBROUTINE F2(T,Y,YP) 

DOUBLE IPRECISION Y(2)) YP G29, 1, DEXP, DSORT, DABS, UA 
DOUBLE PRECISION DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 
DOUBLE IPREGES PONBTET AY CO SDT ERFDVGG 

COMMON /PARAM/DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 
PeUTaY (1) 

V¥EO={V02) 

DTPT=ALPHA*DA1*B1*YCO*DEXP (GAM1*TPT/(1.D0+TPT) ) 
*-ALPHA*TPT 
DYCO=(1.D0-YCO)-DA1*YCO*#DEXP(GAM1*TPT/(1.D0+TPT) ) 
YP(1)=DTPT 


Za) 


100 


206 


YRiC2) =DYCO 

RETURN 

END 

SUBROUTINE FST yy YP) 

DOUBLE PRECISUON (Y(1) PYPIGl) PT SDEXP,DSORT ,DABS,DA 
DOUBLE PRECISION DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 
DOUBLE PRECISION TPT,YCO,ZHC,DTPT,DYCO,DZHC 
COMMON /PARAM/DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 
TPurSay (+) 

YEO=V (2) 

7HCHYy 3) 

DTPT=ALPHAs (DA1*B1*YCO*DEXP(GAM1*TPT/(1.D0O+TPT) )+ 


*DA2*B2* ZHC*DEXP(GAM2*TPT/(1.D0+TPT) )-TPT) 


DYCO=1-D0-YCO#C1. DOTDA 1 ¥DERP(GAMI4TPT/1(1 DOTTPT))} 
DZHC=1.D0-ZHC*(1.D0+DA2*DEXP(GAM24#TPT/(1.D0+TPT) )) 
YP(1)=DTPT 

Ve.G2.)=DYCO 

YPC3)=DZHC 

RETURN 

END 

SUBROUTINE PDERV(T,Y,YP,PW,M,IND) 

DOUBLE PRECISION Y(8,M),YP(M),PW(M,M),T,DEXP,DSQRT 
DOUBLE PRECISION DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 
DOUBLE PRECISION TPT ,YCO, ZHC,DTPT,DYCO,DZHC,DABS,DA 
COMMON /PARAM/ DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 

EP ie Vest 11s) 

VeOanveu 2) 

ZHCsY OP") 

PE CPN SEO A GOLTO47.00 

DTPT=ALPHA* (DA1*B1*YCO*DEXP(GAM1*TPT/(1.DO+TPT) )+ 


*DA2*B2* ZHC*¥DEXP(GAM2*TPT/(1.D0+TPT) )-TPT) 


DYCO=1.D0-YCO*«(1.D0+DA1*DEXP(GAM1*TPT/(TPT+1.D0))) 
DZHC=1.D0-ZHC*(1.D0+DA2*DEXP(GAM2*TPT/(TPT+1.D0))) 
YP(1)=DTPT 

YEPC2Z)=DY CO 

YP (ai) =DZHC 

RETURN 

CONTINUE 
PW(2,1)=-DA1*YCO*GAM1*DEXP(GAM1*TPT/(1.DO+TPT) )/ 


*CORDOFTPT) #*2) 


PW(2,2)=-1.D0-DA1*DEXP(GAM1*TPT/(TPT+1.D0) ) 
PW(2,3)=0.0D0 

PW(3,1)=-DA2* ZHCXGAM2*DEXP(GAM2*TPT/(TPT+1.D0))/ 
CCTPT + H.D0) *#2) 

PW(3,2)=0.0D0 

pw(3,3)=-1.0D0-DA2*DEXP (GAM2*TPT/(TPT+1.D0) ) 
PW(1,1)=ALPHA*(-B1*PW(2,1)-B2*PW(3, le 700) 
PW(1, 2) =ALPHA*DA 1*B1*DEXP (GAM1*TPT/ (1 .DO+TPT)) 
PW(1, 3) <ALPHA*DA2*B2*DEXP (GAM2*TPT/(1.D0+TPT) ) 
RETURN 

END 


SUBROUTINE DFILE(IDCBI,ICRI,IDCBO,NAMEO,ICRO) 
DIMENSION IDCBO( 144) ,NAMEO(3),ISIZE(2),IDCBI(144) 
CALL READF(IDCBI ,IERR,NAMEO) 
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CAEDVCLOSE( IDEBI)) 

PSEZE E48 

Eye r= 

PSECUSd 

CALL OPEN(IDCBO,IERR,NAMEO,0,1,ICRO) 

CALL CLOSE(IDCBO) 

CALL PURGE(IDCBO,IERR,NAMEO, 1,ICRO) 

CALL CREAT(IDCBO,IERR,NAMEO,ISIZE,ITYPE,ISECU,ICRO) 
CALL OPEN(IDCBO,IERR,NAMEO,0,ISECU,ICRO) 

RETURN 

END 

SUBROUTINE RKF(A,Y,N,F,DA,H,HMX,ABSER, RELER, IFLAG) 
*, VARIABLE STEP INTEGRATION OF ODE'S, DTL 11/06/79 


THI Se CODE INTEGRATES. A SYSTEM OF FIRST ORDER ORDINARY 
DIFFERENTIAL EQUATIONS BY RUNGE-KUTTA-FEHLBERG METHOD 
WITH AUTOMATIC ESTIMATION OF LOCAL ERROR AND STEP SIZE 
ADJUSTMENT -- SHAMPINE AND ALLEN (REF 5.16) 


DOUBLE PRECISION Y(N),YTEMP(20) , TEMP(20),R(20) 

DOUBLER PRECISION “K1 (20) /K2(20) /K3.(20) ,K4(20)7K5(20) 
x,K(20),U,RELER,ABSER,HMX,A,B,DA,HMAX,H,HKEEP,ARG, RATIO 
*,DABS ,DMAX1,DMIN1,DSIGN,DSORT,TR,ASV,HMIMA 

REAL DOUT(20) 

COMMON /CT/NPINT,IPINT,NDLAY,IMAX,IBUF(100),IDCBO( 144) 

W=i0De16 

IMAX=0 

TEC RELERY LT.OT0DG ;ORJABSERVLTJ020D0) GO TO 18 

IF (RELER+ABSER.EQ.0.0D0) GO TO 18 

IF(HMX.LE.0.0D0) GO TO 18 

B=A+DA 

ASV=A 

IF (DABS(DA).LE. 13.D0#U*DMAX1(DABS(A) ,DABS(B)))GO TO 19 

HMAX=DMIN1(HMX,DABS (DA) ) 

IF (DABS(H).LE.13.0D0*U*DABS(A)) H=HMAX 

TADIS=0 

H=DSIGN(DMIN1(DABS(H) , HMAX) , DA) 

TR(DABS(B-A) .GT..(1.25D0+107D0#U) #DABS(H)) Go TO 4 

HKEEP=H 

PADI Seu 

H=B-A 

CALE RUAMY | Rap 
' CONTINUE 

Hol cel =n N 

YTEMP(1I)=Y¥(1)+0.25D0*«H*K1(1) 

ARG=A+0.25D0%*H 

CALL F(ARG, YTEMP,K2) 

Dose. =k N 

YTEMP(1)=Y(1)+H*(K1(1)#3.D0/32.D0+K2(1)*(9.D0/32.D0) ) 

ARG=A+H*(3.D0/8.D0) 

CALL F(ARG,YTEMP,K3) 

HOLES TEEN 

YTEMP(1)=Y(1)+H*(K1(1)*(1932.D0/2197.D0)-K2(1) 
**(7200.D0/2197.D0)+K3(1)*(7296.D0/2197.D0) ) 
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ARG=A+H*(12.D0/13.D0) 
CALL F(ARG,YTEMP,K4) 
DOUSHI = wan 

VTEMP GU) =V (1 4H se (Kei Gt?) 
*+K3(1)*(3680. ives pe DO 
ARG=A+H 

CALL F(ARG,YTEMP,K5) 
DOGO = 1 ON 
VIBMPICR)=y (1 4+H*e (-K1(1)4(8.D0/27.D0) +2. D0¥K2(1) 
*-K3(I ere DO0/2565.D0)+K4(1)*(1859.D0/4104.D0) 
=> (is) (11 2D0/40.D0)) 

ARG=A+0.5D0*H 

CALL F(ARG,YTEMP,K6) 

Dee iil = 1> N 
TEMP(I)=K1(1I)*(25.D0/216.D0)+K3(1)*(1408.D0/2565.D0) 
EKA GT ie (2:9 72D0 44104. D0) 0 aA2D04K5 ()) 

ie waren Meat eie 

Pom. ae 

Rit) =K ACI Waren DO-K3(1)*(128.D0/4275.D0)-K4(1) 
Sl esen DO/75240.D0)+K5(1) /50.D0+K6(1)*(2. SNe: DO) 
RATIO=0.0D0 

DOR io R= t 7 N 

TR=DABS(R(1I))/(RELER*DABS (YTEMP(I))+ABSER) 
RATIO=DMAX1(RATIO,TR) 

EPRCURATLO CGT. wieDG) SGONTO S15 

EFA PINT LEB SNDEAN)) GOSTO, 45 

DOS40: E=19N 

BECKeEG! USICS CGT od CD00) eGOeTOmr4 0 

TR CRT) SEOF0 HODGIMGORTO £40 

F((IPINT/IPLOT#IPLOT.EQ.IPINT).AND. 
*UB-A lhl 2D=04)) SIMAX=4 

DE GCLPINT 71 SPEOTS1 PROT. EC SL PINT). AND? 
#(A-ASV.. DT. 0+0001D0)pIMAX=2 
OR A. CNSASMEE TACT OOO0IDO))oy PMAXsa 

TECLMAXMEO A2) ACAULEPOSNT GI DOBO, TERR, - 10) 
IF(IMAX.EQ.2) NPINT=NPINT-1 

DOUT(1)=SNGL(A) 

HMI MA=H*DABS (K1(1)/(K1(1)-K5(1))) 

DOs Sead =n 
DOUT(JJ+1)=SNGL(Y(JJ))+SNGL(HMIMA*K1(JJ)) 

M=N+1 

LE CN SEO. 229 M=aN+2 

CALE LecODE 

WRITE(IBUF,900) (DOUT(JJ) ,JJ=1,M) 
FORMAT(11(1PE13.5)) 

CALL WRITF(IDCBO,1IERR,IBUF,6*M+2) 

NPINT=NPINT+ 1 

GOuTORES 

CONTINUE 

CONTINUE 

DO 514. b=79N 

Y(I)=YTEMP(1) 

A=A+H ‘ 

LP GLADIS. EOmel MCOnTO: 16 


* (4 ae DO/2HG. DOI a8. DOFKR 21) 
)-K4(1I)*(845.D0/4104.D0)) 
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RATIO=DMAX1(RATIO,6.5536D-04) 
RATIO=DMIN1(RATIO, 4096.D0) 
H=0.8D0*H/DSORT(DSORT(RATIO) ) 

PECDABS (Hi) LE el seDO0FUsDABS DAY mIGO TO. 19 
ReCRATT OSE. Wem eee. aren Ss 

PADIS=C 

GO TO '5 

IFLAG=1 

H=HKEEP 

RETURN 

IFLAG=3 

RETURN 

IFLAG=4 

RETURN 

END 

SUBROUTINE SGEAR(N,T,DT,Y,SAVE,H,HMX,EPS,MF, YMAX 
*, ERROR, KFLAG, JSTART, MAXDER, PW, PDERV) 


INTEGRATION OF STIFF AND NON-STIFF SETS OF N ODE'S 
"NUMERICAL INITIAL VALUE PROBLEMS IN ORDINARY 
DIFFERENTIAL EQUATIONS', BY C. W. GEAR (REF 5.15) 


DESCRIPTLONVOF WAR BABLES@IN@ CALL LIST 


Nees ey tO tee Rol ORDERBODE SS 

to. . Lhe INDEPENDENT VARIABLE 

ie COOL aWoLesr) Lom Peel OF LNT EGRAT EPROM ier OemieED T 
(PROBABLY USING SEVERAL STEPS DEPENDING ON H AND HMX) 

Y ... AN 8 BY N ARRAY CONTAINING THE DEPENDENT VARIABLES 


AND THEIR SCALED DERIVATIVES. Y(1,I1) I=1,N CONTAINS THE 
VAGUES SOF y tte SON DYeeney, 1) ES) PROVIDED BY THE 
CALLING PROGRAM ON THE FIRST ENTRY. Y(J+1,1) CONTAINS 
THE J-TH DERIVATIVE SCALED BY H**J/J! 


SAVES pe. WORK AREA OF SIZE 13°°BY°N 
pi. See eos Ai TO" BE AT TEMPTEDOON Tob NEXTS STEERS 
Aw SMALE VALUE OF sh ®SHOULD. BE USED “ON SINE PRS Cai. 
HMAX ... THE MAXIMUM ALLOWABLE STEP SIZE. 
EPS... THE BUCLIDEAN OF THE SINGLE (STREP ERROR GES FiIMAvHsS 


DIVIDED BY YMAX(I) MUST BE LESS THAN EPS. TO ACHIEVE 
THIS, THE ORDER AND/OR THE STEP SIZE IS ALTERED. 
MF ... THE METHOD INDICATOR. 
0 - AN ADAMS PREDICTOR CORRECTOR IS USED. 
1 - A MULTI-STEP METHOD FOR STIFF SYSTEMS IS USED. 
THE JACOBIAN MATRIX MUST BE SUPPLIED BY PDERV 


Pe SAME ASw 41>) EXCEPT. NUMERLGAGS DERE RRANCINGmsU sre 
TO GENERATE THE JACOBIAN. 
VMAX ... VECTOR OF N VALUES WHICH CONTAIN THE MAXIMUM OF 
BAGH! YY tFOUND) SO FAR. NORMALLY SHIT TOU ONSENTRY: 
HRROR@) soeVECTOR OF N= VALUES Oh SINGER, STEP ERRORS OR SYaS 
KFLAG ... A COMPLETION CODE 


ie THES ee OWASeSUCCESEUL 
jo Ter ar AKEN WITH HeHMOIN SUT EPS NG MET 
-2 - THE MAXIMUM ORDER SPECIFIED WAS TOO LARGE. 


OAD A000 10 OOOO OOO O:6-O 2G OO Gy aiere 
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-3 - CORRECTOR CONVERGENGE COULD NOT BE ACHIEVED 
-4 - THE REQUESTED ERROR IS TOO SMALL 
USTAR T=... AN INPUT OINDICATOR 
-| - REPEAT THE LAST STEP WITH A NEW H 
0 - PERFORM THE FIRST STEP (INITIALIZATION) 
+1 - TAKE A NEW STEP CONTINUING FROM THE LAST 
ATahsl Cero TART ELS «SE Tero NO. THE (CURRENT ORDER 
OF THE METHOD. NO IS ALSO THE ORDER OF THE 
MAXIMUM DERIVATIVE AVAILABLE. 
MAXDER ... THE MAXIMUM DERIVATIVE THAT SHOULD BE USED IN 
THE MEHOD. THIS RESTRICTS THE ORDER, AS ORDER 
IS EQUAL TO THE HIGHEST DERIVATIVE USED. 
IF MF=0, MAXDER < 8 IS NECESSARY. 
IF MF=1 OR 2, MAXDER < 7 IS NECESSARY. 
PW ... WORK AREA OF SIZE 2*N**2 
PDERV ... EXTERNAL SUBROUTINE PDERV(T,Y,DY,PW,M,IND) 
THE NAME OF THIS SUBROUTINE WUST BE PLACED IN 
AN EXTERNAL STATEMENT IN THE MAIN PROGRAM. 
IF IND = 0, THEN PDERV IS USED TO CALCULATE THE 
VALUES OF THE DERIVATIVES OF THE DEPENDENT 
VARTABRES (STORED SING Y (1, 1)) l=" N) )sewiTHe RESPECT 
TO THE INDEPENDENT VARIABLE T. THESE DERIVATIVES 
ARE PLACED IN THE ARRAY DY(I). 
Te Mew = sie CHENG ESL NDe =) ile ePDERV STS SUSED TO 
DETERMINE THE ELEMENTS IN THE JACOBIAN MATRIX. 
PW(I,d) SHOULD BE ASSIGNED THE VALUE OF THE 
PARTIAL DERIVATIVE OF THE I-TH EQUATION WITH 
RESPECT TO THE J-TH DEPENDENT VARIABLE. 
PW SHOULD BE DIMENSIONED PW(M,M) OR PW(M,1) 


DOUBLE PRECISION Y(8,1),YMAX(1),SAVE(10,1),ERROR(1) 
*ePERTST(7 , 2-3) ,PWC id, ACS) 
*, BND,D,DBLE,DMIN1,DMAX1,DABS,E,EDWN,ENQ1,ENQ2, 
*, ENQ3,EPS,ERROR,EUP,H,HMAX,HMIN,HNEW,HOLD, 
*, PEPSH,PR1,PR2,PR3,R,R1,RACUM,T, TOLD,HMX,B,TSV,U,DT,S 
REAL DOUT(20) 
COMMON /CT/NPINT,IPINT,NDLAY,IMAX,IBUF(100),IDCBO( 144) 
DaRememREST (2. DOr 4. 5D0 watts 8SD0 10.4 2D0eia D0 
+ eieme OD D0, 2, DO 2 024400 soy 6 SD0g So.3 aD0NW Ue Ucue 
CeCe DUN DO OeD0derlow Dn 2 sol Ups Om CCDU mE UUneu. DUE 
Fete A OUe sy CODO oS oCU0 AFUsCeDU, oun o/DUrae DUR 
Se DOM eNO Oe ODO p01 GG 7D00 Can 3 SD0r 0. U0G2 Gy DUsaeDe, 
mee, DOU |e 0), 2 DUGa DO 0 eo oy D Obey 207 DOs 3 ODO? 
AGT eee cee 0 
U=1.D-16 
IMAX=0 
B=T+DT 
TSV= 
HMAX=DMIN1(HMX,DABS(DT) ) 
HMIN=DMAX 1 (U*T, DMIN1(U,U*HMAX) ) 
IADJIS=0 
H=DMAX1(H,13.0D0*U*DABS(T) ) 
IF (DABS(B-T).GT.(1.D0+10.D0*U) *DABS(H))GO TO 70 
TADUS=1 
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H=B-T 

CONTINUE 

IRET=1 

KFLAG=1 

BE COSTART ILE. 0) GOUTO 140 
DOGO tS -N 

POLO eat, k 
SAVER) = vi Ds) 
HOLD=HNEW 

PAGHVEQ SHOULD) EGOSTC! 430 
RACUM=H /HOLD 

IRET 1=1 

GO TO 750 

NOOLD=NOQ 

TOLD=T 

RACUM=1.D0 

Ee CISTART. Glee )2GO.TO- 250 
GOuTo=170 

BE CISTARTOES S317 )0GO TO 160 
NQ= 1 

N3=N 

N1i=N*10 

N2=N1+1 

N4=NxX*2 

N5=N1+N 

N6=N5+1 

N= 1 Z25N 

DO a5 eol=N 
SAVE(N7 +00, 1)=0.0D0 
CALE PDERVCE,Y,SAVE(N2,1),PW,N3,0) 
DO, 1500L=47N 
M2) =SAVE (NGI 7 1) 4H 
HNEW=H 

K=2 

COBTOe G00 
IF(NQ.EQ.NQOLD) JSTART=1 
TeTOnD 

NO=NQOLD 

K=NQ+1 

Gator I20 

IF(MF.EQ.0) GO TO 180 

TE (NO.GTV6)- GO TO 190 


WiGOATOS C22 Ie oo eee 245 22 oe Zoe NO 


BE CNO SGT 87) GO TOM 90 


GORTOU? bie 2 eee 145 245 ore eNO 


KFLAG=-2 

RETURN 
Ruly=—1eD0 

GORTOR2 30 
A(1)=-1.D0/2.D0 
Als) =—1)2D0 72 aD0 
GORTOSZ30 
A(1)=-5.D0/12.D0 
A(3)=-3.D0/4.D0 
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ACA) == 1. D0/6-.D0 
GO TO 230 

214 R= =8"po 738 -D0 
Pee = ee Oa De 


Ae) == eDOY Sept 
AM 5)Ss1.D0/24 530 
GORTORZS0 


215 Pu == 25 120 Ie Or De 
Aes) =—25" D072 245 Do 
Rueay=—3 5eD0/7 2200 
AUSy==5 D0 /483D0 
PAG) == SODA 720R 0DO 
GOnTOG2 30 

216 Ati =—955DO728 82 De 
Pars) 3 7 20641 20) DO 


= 

ACS) ==5 2D0/8.D0 
MGS) =at22D0796 2D0 
Gia |) D0740 D0 
Rey == 15 D092 07.00 
GO?TO7230 


2457 A(1)=-19087.D0/60480.D0 
A(3)=-49.D0/40.D0 
Pet) =—=2080D07 270 .D0 
A(5)=-49.D0/192.D0 
AC6)}=—750D07 144 (D0 
ACT) =-7. DOs 14407D0 
A(8)=-1.D0/5040.D0 
GO8TO +230 

Dior PCiss DO 
GO TO 236 

222 Bu) =s20D0 7/3. D0 
MGs j== F200 739D0 
GOnTOR230 

228 Me=-6 5007 Ir. DU 
AOS) ==67D0/7 19eD0 
Aiay=—-49D0/11.D0 
GO-TO8230 

224 eM =— 12 D007 25.00 
AGs == 7 eD0 £10. D0 
A(4)=-12D0/5<D0 
A(S5)=-1.D0/50.bD0 
G@ATO 72350 

Dole N= he ODO 27 ep 0 
A(3)=-225.D0/274.D0 
A(4)=-85.D0/274.D0 
A(5)=-15.D0/274.D0 
A(6)=-1.D0/274.D0 


COBO eZ 20 

226 A(1)=-180.D0/441.D0 
A(3)=-58.D0/63.D0 
A(4)=-15.D0/36.D0 
A(5)=-25.D0/252.D0 
ACG) ==37D07252 200 
B07.) =) 2 D0 /s7 64. D0 
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260 


270 


280 
290 
300 


ord 
320 


330 
340 


S530 
360 


S70 
380 


K=NO+1 

I DOUB=K 

MTYP=(4-MF) /2 
ENQ1=0.5D0/DBLE(FLOAT(NQ) ) 
ENQ2=0.5D0/DBLE(FLOAT(NQ+1)) 
ENQ3=0.5D0/DBLE(FLOAT(NQ+2) ) 
PEPSH=EPS 
EUP=(PERTST(NQ,MTYP, 2) *PEPSH) **2 
E=(PERTST(NQ,MTYP,1)*PEPSH) **2 
EDWN= (PERTST(NQ,MTYP, 3) *PEPSH) **2 
IF(EDWN.EQ.0.0D0) GO TO 780 
BND=EPS*ENQ3/DBLE(FLOAT(N) ) 

I WEVAL=MF 

GOPrOe250 680) 1 RET 

T=T+H 

DOS260 I= 27K 

DGS2602 01-0 eh 

J2=K-Ji+J-1 

DOEZ60 e1=1,N 

Vito waa v. Ils oey Cio 1.1) 

Den 270 lea. N 

ERROR(I)=0.0D0 

Do, 230 b=1,3 

CAlMeEDERV (ly SAVE(N2 “1)PWON3, 0) 
PEAT WEVALOISE 1b IGO. TO 350 

MEME ZEOs2) GO TOr3 10 

CALL PDERV(T,Y,SAVE(N2,1),PW,N3,1) 
R=A Gi) +H 

DOwZo0) t=1,N4 

PW(I)=PW(I)*R 

DO 300 W=1 ,N 
PW(1*(N3+1)-N3)=1.D0+PW(1*(N3+1)-N3) 
IWEVAL=- 1 

CALL MXINV(PW,N,N3,24N3,J1) 

Pet eGra0 )eGO ITO S350 

GO TO 440 

BO 320 l=1,N 

SAVE GS ar) =v (1 wh) 

BOBSA0 7I=1.N 
R=EPS*DMAX1(EPS,DABS(SAVE(9,J))) 
Vien Wied oR 

D=A(1)*H/R 


TORU EPDERVCUT AY, SAVE NGaIy , PW Nas 0) 


DOr 360. b= uN 
PW(I+(J-1)*N3)=(SAVE(N5+I,1)-SAVE(N1+1I,1))*D 
Vit Sa) =SAVE(S 7d) 

GOLTOe290 

IF(MF.NE.0) GO TO 370 

DOSs6ORL dN 

GAVEGS 1h) Vice si) SAVE Na ieg e 
GOErOeaI0 

DO 380) 1=15N 
SAVE(N5+1,1)=Y¥(2,1)-SAVE(N1+1,1)*H 
DOELUO0 sieriN 


PS . Te a 
=i} RAN ; 


: : 
a : 


390 


410 


420 
430 


440 


460 
470 


490 


500 


BS hU 
520 


530 


540 


55.0 


1) *N3) *SAVE(N5+J3,1) 


N 

1) +A( 1) *SAVE(9 1) 

I)-SAVE(9,I) 

ROR(1I)+SAVE(9,1) 
E(9,1)).LE.(BND*YMAX(1))) NT=NT-1 
CONTINUE 

LE(NTSLESO} GO TO 290 

CONTINUE 

IADJS=0 

T=TOLD 

IF((H.LE. (HMIN*1.00001D0)).AND. 

*( CIWEVAL-MTYP) .LTU-1))} GO TO 460 

IF((MF.EQ.0).OR. (IWEVAL.NE.0O)) RACUM=RACUM#0.25D0 
IWEVAL=MF 

IRET1=2 

GOnTOUT SO 

KFLAG=-3 

PORES ORI = TaN 

DOL400 30= 15K 

WCU) =SAVE (I, 1) 

H=HOLD 
NQ=NQOLD 
JSTART=NO 
RETURN 

Bs 0) 6DG 
DOS BiaitEN 
D=D+ (ERROR(I) /YMAX(I))**2 
O TO 540 

Omron 520 


Vi@jmm)=V CJ, 1) +A( J) ERROR GI» 
KFLAG=1 

HNEW=H 

HeeipOUs bE: 1) GO-TO 7558 


I DOUB=IDOUB- 1 


SPR CMDOUB.GT.1t)) GOO TO@7/00 


DO Sa0 t=a1 N 

SAVE(10,1)=ERROR(1) 

GOn tO 00 

KFLAG=KFLAG-2 

TADIS=0 

IF(H.LE. (HMIN*1.00001D0)) GO TO 740 
T=TOLD 

LTE CKRERAG TEE. —5)"GO TO 720 

PR2=(D/E) **ENQ2*1.2D0 

PR3=1.D+20 
IF((NQ.GE.MAXDER) .OR. (KFLAG.LE.-1)) GO TO 570 
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D=0.0D0 
DOC O rad eN 

560 D=D+( (ERROR(1)-SAVE(10,1)) /YMAX(I))*«2 
PR3=(D/EUP) **ENQ3*1.4D0 

570 PR1=1.D+20 
DR INCw EEN) GOmTO 6591) 
D=0.0D0 
DOSSE 0 era iN 

580 B= Dewy kh ale) A YMAX Gls) ko 
PR1=(D/EDWN ) **ENQ1*1.3D0 

590 CONTINUE 
PRCPRZ bE. PR3) GO TO 650 
IF(PR3.LT.PR1) GO TO 660 

600 R=1.D0/DMAX1(PR1,1.D-04) 
NEWQ=NO- 1 

610 I DOUB=10 
PUmGnrnAG- KOA AND a(R. bis (15 1.D0)))).GO TO" 700 
IF (NEWQ.LE.NQ) GO TO 630 
BORG? 0iI= 1 N 

620 Y (NEWO+1,I)=ERROR(I)*A(K) /DBLE(FLOAT(K) ) 

630 K=NEWQ? 1 
Hathe LAG. EO. i GOTO (670 
RACUM=RACUM#R 
IRET1=3 
GO TO 750 

640 IF (NEWO.EQ.NQ) GO TO 250 
NO=NEWO 
GOsTO 70 

650 IF(PR2.GT.PR1) GO TO 600 
NEWQ=NQ 
R= 17D0/MMAX1(PR2, 1.D-04) 
GO ETO =6710 

660 R=1.D0/DMAX1(PR3,1.D-04) 
NEWO=NO+1 
GO TO 610 

670 IRET=2 
R=DMIN1(R,HMAX/DABS (H) ) 
H=H*R 
HNEW=H 
IF(NQ.EQ.NEWQ) GO TO 680 
NO=NEWO 
GOuroni 70 

680 ~~ Ri=1.D0 
DOmoS0) d=2> Kk 
R1=R1*R 
DOeGSOMI =e N 

690 WiGuele= Yue) ER | 
I DOUB=K 

700 Oe OMe N 

710 YMAX (I )=DMAX1(YMAX(I),DABS(Y(1,I))) 
JSTART=NQ 
Pri PINTeEEYNDLAY,)) GO) TO 7 14 
DOM Tie Juan : 
TPUUSAVEUNG Ud Bl) *Y (2,00 )),Gl.0.0D0).. OR. 
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Peon] aI te sO. OD0))) GOuTO, 712 
TRCOLPINT/IPLOT#IPLOTSEO.IPINT).AND. 
FRC eject Or, OOOO) eel MA x= 1 
IFC GLRINT/IPLOT41I PLOT «EQ. 1PINT) SAND. 
HEGT=TSV)°bT.020001D0).)) TMAX=2 
IF (IMAX.EQ.2) CALL POSNT(IDCBO,1IERR,-1,0) 
IF(IMAX.EQ.2) NPINT=NPINT- 1 
BOUT (1)=SNGL(T) 
DOV713 KK=1,N 
P13 DOUT(KK+1)=0.0 
SaeDERS(Y (2 oy odd SAVE (NZ todd) )) 
NQQ1=NO+1 
DOS SeKK=12N 
DOM ae = 1" NOOK 
711 DOUT(KK+1)=DOUT(KK+1)+SNGL(Y(3,KK)*S*#*(J-1)) 


€ DOUT (RK) =SNGDCY (1) KR) +S ehey G2) KK) 
715 CONTINUE 

M=N¢+ 1 

IF(N.EQ.2) M=N+2 

CALL CODE 


WRITE(IBUF,900) (DOUT(KK) ,KK=1,M) 
900 FORMAT (114(1PE13 75) ) 
CALL WRITF(IDCBO,IERR,IBUF,6*M+2) 
NPINT=NPINT+1 
GOeTrO 3714 
P12 CONTINUE 
mee DO At oad Sanan 
716 SAVE(NG a yl = vee!) 
i GLAbdS BO. OGG TOus0 
RETURN 
720 PeON@ sho.) GOmTO 7380 
CALL PDERV(T,Y,SAVE(N2,1),PW,N3,0) 
R=H/HOLD 
DORIC Onl=1.N 
VY (erie = SAVE, Ey) 
SAVE (Qal V=HOED*SAVE(N 1 +1 a) 
Tea8) ViG2 I= SAVE OWL JER 
NQ=1 
KFLAG=1 
GOsTOm LO 
740 KFLAG=—- 1 
HNEW=H 
' JSTART=NQ 
RETURN 
750 RACUM=DMAX1 (DABS (HMIN/HOLD) , RACUM) 
RACUM=DMIN1(RACUM, DABS (HMAX/HOLD) ) 
R1=1.D0 
DON A600 Sd =—2y k 
R1=R1*RACUM 
DOM 60nL— iN 
760 ¥ (3, DJESAVE CI, 1) #R4 
H=HOLD*RACUM 
DOM Or =7N 
a ViGl P= SAVE G1) 


TUGs 


hi : 
a 
oe | =. agate cat tae Th yee 
iy i we hd iy = > ii yi woe 
fa § ’ i, : 7 ; 
| (tc wr P vim 


a) 
_ ; wN. 
iu 


I DOUB=K 
GOMTOm (130 R250°640)) TRET4 
780 KFLAG=-4 


IF (DABS(B-T) /GT.(1.D0+10.D0*U)*DABS(H)) GO TO 470 


T=B 

GOV TO: FIG 
(& GO TO 470 

END 
Cy DECK MZINV 
C 


C 
c SUBROUTINE MXINV 
G 
C PURPOSE 
C INVERSION OF A REAL SQUARE MATRIX. 
c 
C USAGE 
fe CALL MXINV(A,N,M,M2,IER) 
G 
C DESCRIPTION OF PARAMETERS 
CuK MATRIX TO BE INVERTED. 
Gan - A MATRIX DIMENSION (UPPER LEFT CORNER) 
CM ~ ROW DIMENSION OF A IN CALLING PROGRAM 
C M2 - 2%M # OF COLUMNS IN A MATRIX IN CALLING PROGRAM 
C.RER - ERROR FLAG. 
G =1, NORMAL RETURN. 
C =-1, MATRIX A IS ALGORITHMICALLY SINGULAR. 
‘s 
c MATRIX A IS DESTROYED, AND IS REPLACED BY ITS INVERSE 
G 
G SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
G DABS (FORTRAN ) 
© 
C REFERENCE 
‘e MATRIXALGEBRA OG ALGEBRAISKE LIGNINGER. 
C VILLADSEN,J. 
G 1968, POLYTEKNISK FORLAG, COPENHAGEN, DENMARK. 
c 
C 
SUBROUTINE MXINV(A, N, M, M2, IER) 
DOUBLE PRECISION A(M,M2), DET, OPSTO, DABS 
INTEGERS gal RR eo oh ee NINE een Te? 
Ni Nek 2 
DET = Ww, D0 
NE le aNe al 
DOs 4.0 ial aeleaeN 
DO 40 J = NPI, NT2 
Aleem ODO 
[ee Giese eN) . 407 230) Be0 
30 AG? dee ae, 0 


40 CONTINUE 
pO 140 f= 15 N 


? 
co 
J 


} 7 ‘ q ae 
> ee a z se. oe 


! u 
: neo, 


nd ae ata Fe 


only 


50 


60 


PER CDARS CAG meer OCS TON 6Or G0r 


ra 


POR me S 


OPSTC 


LP i= T +4 
OPSTO\ = DABS(A(1,1)) 
Keet 
Donec Ot eea ees 
ce al 
OPSTO1= DABS (A(G 1 )})} 
CONTINUE 
eK = a7 Oe 90 
DET = -DET 
HOME 0M  =ehe NT? 
OPSTOR=TA (le) 
Phe eye AK 
ACK) ="OPSTO 
Poe Atk) e100 
OPSTOM=Fie DOV FACT SL) 
DET = DET / OPSTO 
DOetiONS. =als NTZ 
Riles). =A ope 
DOOR TESS EN 
Race ae. 20, 
OPSTOe=8A CIs) 
IOP GhS1e: ie coe Eee. 
ACK) GHA Gee) 
CONTINUE 
DO Omia= a7, aN 
Hom oO Sree N 
NI =N+J 
ee seACLENG) 
CONTINUE 
IER=1 
RETURN 
IER = -1 
RETURN 


END 


50 


les) (COC) (BONO) 1) 


TN 


800 
805 


801 


704 


FO5 


706 


SiH 


PS) 


THIS PROGRAM IS USED TO CALCULATE THE AUTOCORRELATION 
COEFFICIENTS FOR A SERIES OF TEMPERATURE MAXIMA. 

THE AUTOCORRELATION COEFFICIENTS CAN BE CALCULATED 
BASED ON EITHER THE AMPLITUDE OR THE FREQUENCY OF 

THE TEMPERATURE MAXIMA, 


Arex 
PROGRAM AUTOC(,30500), CALCULATION OF AUTOCORRELATIONS 
DIMENSTONGIDCBI (144) NAMED (3)))T( 10000) ,AC(50) > 2(4) 
DIMENSION IBUF(100),NAMEO(3),IDCBO( 144) ,2ZS(4) 
DOUBLE PRECISION T1BAR,TKBAR,VART1,VARTK,VAR 
x ,DBLE,AUTCV 
LUNW=LOGLU(ISIS) 
LUNR=LUNW 
CALL IFMTE(LUNW) 
CALL IEROE(LUNW) 
WRITE (LUNW, 800) 
FORMAT(" ENTER INPUT DATA FILE NAME, AND ICRI") 
READ(LUNW,805) NAMEI,ICRI 
BORMAT GAZ aries) 
WRITE(LUNW, 801) 
FORMAT(" ENTER OUTPUT FILE NAME, AND ICRO") 
READ(LUNW,805) NAMEO,ICRO 
CALL (DFILE (1 DCBO /NAMEO, ECRO) 
CALL OPEN(IDCBI,IERR,NAMEI,0,1,ICRI) 
CALL READF(IDCBI,IERR,IBUF,4) 
GALE CODE 
READ(IBUF,*) NINT 
WRITE(LUNW,704) NINT 
FORMAT(" NINT = ",15) 
WRITE (LUNW, 705) 
FORMAT(" ENTER NINT, AND K (# OF AUTOCORRELATIONS) ") 
READ(LUNW,*) NINT,K 
WRITE(LUNW, 706) 


FORMAT("ENTER TYPE OF A-C, 1=TIME(FREQ), 2=TEMP(AMP)") 
READ(LUNR,*) ITYPE 

KP1=K+1 

GALE CODE 


WRITE(IBUF,910). KP1 
CALL WRITF(IDCBO,IERR,IBUF,3) 
FORMAT (16) 

’ CALL READF(IDCBD, TERR, BUR, 26) 
GALLILWRITE (PDCBO/LERR, IBURy 26) 
CALL READF(IDCBI,IERR,IBUF,20) 
CALL WRITF(IDCBO,IERR,IBUF,20) 
GALL READE (IDCBI, TERR; IT BUER39) 
GAUG WRITE LIDCBO, IL ERR, DBUR 73S) 
DRCLTYPESEO; 2) GO TO 1¢@ 

DOneSe rl =feNINT 

CAMESREADE (LUCBI (TERRE LSUR, 26) 
CALESCODE 

READGIBUR, m2 (01), il=174) 
Pe) = 22) Cale) 


an . ¥ 
Ree +. was bad 


« 4 nase aaa c* 


BY os eee C5K2 
ages f eas) SGb2 <lLh 
7 bye ad! 


aa | =, 4 
7 > 
# 


1 
20 


30 


40 


50 


60 


920 


80 
g00 
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GOSTOR ZO 

DORIS CT =1eNINT 

CALL READF(IDCBI,IERR,IBUF,26) 

CAEL CODE 

READ CUBUP P+) a GZiCU Tent t=) 

Oe =z (3) 

CALL CLOSE.CIDCBI ) 

DOme On hayek 

NSUM=NINT-IK 

T1BAR=0.0D0 

TKBAR=0.0D0 

DO 30 ISUM=1,NSUM 
T1BAR=T1BAR+DBLE(T(ISUM) ) /DBLE(FLOAT(NSUM) ) 
TKBAR=TKBAR+DBLE(T(ISUM+IK) ) /DBLE (FLOAT (NSUM) ) 
VART t=0.0D0 

VARTK=0.0DO 

DO 40 ISUM=1,NSUM 

VART 1=VART 1+ ((DBLE(T(1SUM) )-TIBAR)#*2)/ 
*DBLE(FLOAT(NSUM) ) 
VARTK=VARTK+ ( (DBLE(T(ISUM+IK) )-TKBAR) **2) / 
*DBLE (FLOAT (NSUM) ) 
VAR=(VART1+VARTK) /2.0D0 

AUTCV=0.0D0 

DO 50 ISUM=1,NSUM 

AUTCV=AUTCV+ (DEGECTCISUM) /=TIBAR) * (DELE(T(1SUM+IK) ) 
*-TKBAR) /DBLE(FLOAT(NSUM) ) 
AC(IK)=SNGL(AUTCV/VAR) 

ACO=1. 

TO=0 

CALL CODE 

WRITE (LBUF,920) 10,ACO 

CALL WRITF(IDCBO,IERR,IBUF,9) 
WRITE(LUNW,920) I0,ACO 
PORMAT (14 1X (PE 13.5) 

DOmo Omelet ik 

GALL CODE 

WRITE(IBUF,920) IK,AC(IK) 

CALL WRITF(IDCBO,IERR,IBUF,9) 

WRETE (LUNW, 920) TK, ACUI) 

FORMAT (6(1PE13.5)) 

GALI. CODE 

WRITE(IBUF,900) T1IBAR,TKBAR,VART1,VARTK 
CALL AWRETE (CIDCBO STDERR, 1 SUP) 26, 
WRITE(LUNW,900) T1BAR,TKBAR, VART1,VARTK 
CALL LOCE CIDGBO, TERR, 1 REC, IRE TOrP SHC) 
ITRUN=(JSEC/2)-(IRB+1) 

CALL CLOSE CILDGBO, TERR,ITRUN) 

STOP 

END 

SUBROUTINE DFILE(IDCBO,NAMEO,ICRO) 
DIMENSION IDCBO(144),NAMEO(3),ISIZE(2) 
TSE.ZE GD a 

ITYPE=3 

ESECUET 


CALL OPEN(IDCBO,IERR,NAMEO,0,1,ICRO) 

CALE CLOSE. CIDEBO? 

CALL PURGE(IDCBO,1ERR,NAMEO,1,ICRO) 

CALL CREAT(IDCBO,/1TERR,NAMEO| [SIZE ITYPE,I1SECU,TCRO) 
CALL OPEN(IDCBO,IERR,NAMEO,0,ISECU,ICRO) 

RETURN 

END 
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20 
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PROGRAM USED FOR MULTIPLICITY AND STABILITY ANALYSIS 


FOR 3 ODE SYSTEM FOR BETA2 VARYING WITH OTHER PARAMETERS 
HELD CONSTANT. 


EXTERNAL FT 
DIMENSION EVAL 
COMPLEX*16 W(3 


EQUIVALENCE 


Ro SOE sony SS ON ESS CO, eso) 
) 


(Wit) SEVAGT 1. 1)) 


COMMON /BLK1/ DA1,DA2,B1,B2,GAM1,GAM2,ALPHA 


READ ( 


BCU) abs 


FORMAT(6D10.5) 
5,905) NDATA 
FORMAT(I2) 


READ ( 


INR=7 


TeDAZMBANGAMiIPGAM2 ALPHA 


MM=NDATA- (NDATA/7 ) #7 
NOALL=NDATA/7+1 


IF (NDATA.EQ. 


(NDATA/7)*7) NOALL=NDATA/7 


DO 200 IOALL=1,NOALL 


TH@GROALESEOUNOALL) .AND ACNDATA.NE. UNDATA/7) #7 


ees 805) 
FORMAT ( ry TLD) 

WRITE(6 810) DA1,DA2,B1,ALPHA,GAM1,GAM2 
FORMAT (T20,'DA 


eEZOF 
x EZOU 
DOO 
READ ( 
DO 20 
TSS (I 
EPS=1 
EPS2= 


ETA=0. 


NSIG= 
N=5 


"GAM1 = 


Bl = 


0 IDAT=1 
5,900) B2 
TJ=1,5 


J)=((B1+B2)/5.0D0)*DFLOAT (I 


> OBevldls 
1.0D-04 
10D0 

2 


ITMAX=50 


IER=0 


1 = pplPOl eet 1, BDAZF ~= oie, Dial 


aD Preto SU ALPHA=" "9 D11.454, 
Dr eter SH UGAMIy =e D1 104) ////) 
,INR 


CALL UERSET(2,LEVOLD) 
CALL ZREAL1(FT,EPS,EPS2,ETA,NSIG,N,TSS,ITMAX, IER) 


NSS=0 
DO 30 


VER UGLSS Glu eebe 


IJ=1,5 


NSS=NSS+1 


TSS(NSS) 


SESS Gh) 


FAC 1=DEXP(GAM1*TSS(NSS)/(1.D0+TSS(NSS) ) ) 


)) 


INR=MM 


ee 


= eB i+ B2)-*0.,.0 1 D0 


C20)}sORTGTSS URI SGT A(R Boe): GOsrOsst 


e'Y =',D11.4 nen 


FAC2=DEXP(GAM2*TSS(NSS)/(1.D0+TSS(NSS) ) ) 


SNS ee Oe ee 
ZSS(NSS)= fas .DO+DA2*FAC2) 


Guugu ee er BRUSS(NSS)) VSS(NGSimZSSUNGoIy 
WRITE(6,820) b2,TSS(NSS) ,¥SS(NSS) , 25S (NSS) 
FORMAT (' ',T10,'B2 =| Digi AT oC eas eee 6, 


2 =) bi :8) 
BB(2)=BB(2) +BB ( yas 


ves a PS eee ek Waa 


‘OLITEATS Cid PRIS ee 


eee . 
Manet omives aT Sati isa’ SRD, Ae ee 
a : Rica, (Ope sy ~ PET She ny dave wore! 
Oe: see 
‘Tt. tO GaVRL SI i) open ng 
ch 18. Sag fe Vea 
MAS ta, SABE ef ibs Bee 
is ‘Bi ott el 
ae “14 (afi ie Pres 
- FEST AMACS 
Your el” 
% hit at Aa 
‘ht 4 Tete (TAS 
ete a i, Da, APACE RI 
Ac AT. OO Oo 
| Wid fad Lie yO AGS Pee 
Leos ya AT ae 
. o, 0? Se emagoe 
P s 1) S od FAT tA 
beret : dc? .O0T aeos 
Lie Mai the. ,<OaTe: 
ss fdo") & FeAD” OST 
Cui, GAT 20} aa 
.£e (600,24 tAsn es 
o hed, Glia 
, ; inf oo €v< See) Ge leliriges 


aged, TeB5@ 

he-@0 7 *taoa 
GO07 . (ATE 

‘ 2e3igH 
: ao tae 
NS wAAMDT 

iwegy 

*J0veu, §)TSaae aS 
242 0) ( Ame aad... 
O=2eH . 


Pom Tweety: eT 


Of oF Goi fi mhi Ww). LO, (bier) mo 


mae \ 


2) “a2 QE Cuts’ 


_ 7 
20%: 
on 
ve 
_ . 
7 AS 


sing ie _ 


nea * rig aa 


Ries nd ¥ oe 


S20 


810 


Aas 


~— 


WRITE(6,825) (BBC L)e 
BRORMAT (110, “TRAGHS=irs 
*D11.4,T61,'DET =! ie 
WRITE(6,830) (EVAL (1, 
ORIEN We, HR GE ie 

*T60),D11.4:) 
WRITE(G,835)) (BVAL (2, mA ne 4°°3)) 
FORMAT(T23,'IMAG' 730, ie cee Dial al GON dies tee, 7a) 
CONTINUE 
CONTINUE 
CONTINUE 
STOP 
END 
DOUBLE PRECISION FUNCTION FT(T) 

IMPLICIT REAL*8(A-H,0-2Z) 

COMMON /BLK1/ DA1,DA2,B1,B2,GAM1,GAM2 
T2=T 

eT eae Os D Ona = 0 Daa 0) 
F(T.GE.B1+B2) eae ant aa 
T3=DEXP(-GAM1*T2/(1.0D0+T2) ) 

T4=DEXP(- ee .OD0+T2) ) 
FT=DA1*B1/(T3+DA1)+DA2*B2/(T4+DA2)-T2 
RETURN 

END 

SUBROUTINE EVALUE(W,BB,T,Y,2Z) 
IMPLICIT REAL#8(A-H,O-Z) 

DEMENSTONGR (3c 05 WKS) BB CS) 

COMPLEX * 16 woae 22(3,3) 

COMMON /BLK1/DA1, DA2,B1,B2,GAM1,GAM2,ALPHA 
IER=0 


$ \PUNe 26 ke SOMES DE TS ATR =, 
) 


eater 
REAL VSE3Z0N1PDId 40045 D1 ond), 


RJ(1,1)=-1.D0-DA1#DEXP(GAM1*T/(T+1.D0)) 

RevGtee yy =0" 0D0 

RJ(1,3)=-DA1#DEXP(GAM1*T/(T+1.D0) )*Y*GAM1/(T+1.D0)*«*2 
Regen) 02 000 
Ro(2,2)=-1.0D0-DA2*DEXP(GAM2*T/ (T+1.D0)) 

Ro (2 53) =-DAZ*DERPAGAM2Z47/ (T+ 0. D0) ) #2GAM2/4 7 2 DO 32 
RJ(3,1)=ALPHA*DA1*B1*DEXP(GAM1#T/(T+1.D0) ) 
RJ(3,2)=ALPHA*DA2*B2*DEXP (GAM24«T/(T+1. mee 


GG aA LPHA® U= Bd RUNG) Bose se OD 

ys eal PRE Ro OREO 

SMIENOR=RG( 1) PY eR (22) RON 2 eG G2 Jee Rg ees Ree 
¥“RI(3, 2) RO (273) FRU 1) eR sy SER U3 yl eRO i 3) 


SOET=RI GP) 4R I (242) + Roa Re ore Rn gers eRe 


PERC 2) RIA, OARS) deORD Wepe) eRe tbe cos eae 
PERI 3) eRe Gs, 2) ROC) RO ee eRe Os Rose 2a) 
BB(1)=TRACE 
BB(2)=-1.D0*SMINOR 
B(3) =DET 
GALL -EIGRE(RG 35.0) Gy W, ooo pn, LER 
IF(IER.LT.32) RETURN 
WRITE(6,810) IER 
FORMAT(//,1T20,'IER FROM EIGRF = ',14) 
RETURN 
END 
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Appendix C Rotameter and Flowmeter Calibrations 


SHOwWhmaiMperTgures CG. 1 )to C.>oeare calibration curves for the 
two rotameters and the two linear mass flowmeters. The 
rotameters and mass flowmeters were calibrated at both 20 
psig and at 30 psig. However, unlike the rotameters, the 
mass flowmeters are not affected by changes in pressure. 
Thus the calibration curves in Figures C.4 to C.5 are valid 
over a wide range of pressures (1 psia to 250 psia). In all 
cases the flow rates are given at the standard conditions of 


zero degrees Celsius and one atmosphere. 


Flow Rate, em ”/sec (SER 


Figure C.1 
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Rotameter Calibration Curves for Oxygen 


steel Ball 30 psig 
Steel Ball 20 psig 
Glass Ball 30 psig 
Glass Ball 20 psig 


30 60 90 120 
Rotameter Reading 


(Matheson #600 tube) 
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150 


= ae et 
4A He 30 psig . 
M@ He 20 psig 


Flow Rate, om ”/sec (STP) 


20 AO’ 60. eO HOO 
Rotameter Reading 


Figure C.2 Rotameter Calibration Curves for Helium 
(Matheson #610 tube) 
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® Ne 30 psig 
* Ne 20 psig 


wv) es) aN 


Flow Rate, em ”/sec (Si 
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Figure C.3 Rotameter Calibration Curves for Nitrogen 
(Matheson #610 tube) 
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Figure C.4 Calibration Curves for Matheson Model ALL-100 
Linear Mass Flowmeter (0-100 SCCM) 
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Figure C.5 Calibration Curves for Matheson Model 8116-0153 


Linear Mass Flowmeter (0-5000 SCCM) 
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Appendix D Gas Chromatograph Calibration 


PieravheceD | touD Searergiven thescalibration datawtor the 
gas chromatograph, in which the response factors of carbon 
dioxide, oxygen, and carbon monoxide relative to nitrogen 
are presented. No calibrations were performed for water, as 
it was not present in any of the reaction mixtures. These 
response factors are necessary in order to convert raw data 
inches form of -areaspercent(orsin arbitrary units) Inte 
actual compositions in the form of mole fractions. The 
calculation procedure which must be followed in order to 


Becrormuch i Ssyconversion 1s sguvenetollowing Table D3. 


Hable D.1" GG Calibration Data for CO,-N-] Mixtures 


GOP Ac eula | CO,, Area Rees COD NepAcE Ua / 
N, Ne COSTS Area 


ORO 7 O02 OLS 4.5 0.856 
Oe ah tele 0.864 
ORS 77 0.860 
0.01614 0. On 8 w2 0.862 
0 01870 0, 8.63 


OOO 0.864 
0.01883 Ore 50 0.864 
OO 21:65 OFO70 


0.02584 0.860 
OmO2 221 0.02566 0.866 
O02 56 02869 
0.02547 ONG 2 


O00 tae 
Oi et 
0.03104 


0.04357 0.865 
0S037-68 0.04343 0.868 
0.04341 0.868 
0204325 Cele 
U0 65,95 0. 8%3 
5572 0.0637] O28HS 
0.06342 0.880 


Average Response Factor O.867 
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Table D.2)'GO Calupration Data for O7-N. Mixtures 


O,, Actual O,, Area On /NaeAGeia) 
N, N, Oe aE Area 


~UOG 
0.2000 Oe 
1.065 


ie SENS) 
IOS 
069 


Ors S21 


On 324 


Average Response Factor 


Tables .3 °tGG Calibration Data for CO-N, Mixtures 


co, Actual oh Area RF = CO/N, Actual 
COCO Aa Area 
0.01866 62-04 52a Bee 
0.01491 e252 


0.01808 1.243 
0.02248 OF Gi eug e248 
O72 Oe S tegtey 


0. 030u-9 LPararas) 
0.03685 OU SOR ieee 
OG Snes ier 23 


0.07638 Lege 

e026 9 OF C/6a).6 13220 
O°. 0 7-57 4 i226 

ene es S 

Urner 396 ea oe 1.246 
1589] oe 
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GC Calculation Procedure 
In order to perform the calculations it is necessary to 
assumeguhatethetthermal@conductivity of a mixturesofrO, , Ne 
and CO is a linear function of the pure component thermal 
conductivities and mole fractions. Using this assumption, 
the split ratio (SR) between the two parallel branches of 
the column arrangement is give by: 
Sheoenmea (OQ, ) ts AteaiN,)) tuAreaCO) (irom oA MiSs) (Doe 
~ Area(O,+N,+CO) (composite peak from P-S) 
where M.S. stands for molecular sieve, and P-S stands for 
Porapak T - Spherocarb. Once the split ratio has been 
calculated, the amount of CO, which has been held up in the 
5A molecular sieve column can be determined from knowing the 
amount of CO, which passed through the Porapak T - 


SoResocathD pat, ol. the covumn, 2.e.:: 
Area\ GO 7@Erome5A M lS?) "=(SR*¥Area (CGO, from, P-S) (De2) 


Ties calculation of the mole fraction of each component can 


now be performed by the following procedure: 


CD = Area(O,)*RF(O,) + Area(N,)#RF(N,) + Area(CO)*#RF(CO) 


+ SR*Area(CO, from P-S)*RF(CO,) (Des) 
MolémFraction O9e= Atea (Or )¢RE (07) 6D (Ded) 
Mole Fraction N, = Area(N,)*RF(N,)/CD (D.5) 
Mole Fraction CO = Area(CO)*RF(CO)/CD (D.6) 


Mole Fraction CO, = Area(CO, from P-S)*RF(CO,)*SR/CD (D.7) 
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Where Rr Stands) for response factor wivAsnumerical example of 
this calculation procedure is given below, where it is 
Scie imeNatwas muxt Ure sot eco s COS OlanGuNoinhas wcauseartne 


pouegrator £O produce the. following raw data: 


Component Area% RF(irom.Tables#D.1-D. 3) 
OFA N2 =CO a5 — 

Ges 10 O28o7 

0; 20 1,067 

N, 30 ihe le 

CO 5 hos Poche 


Brome EQuatwon Dl theesplit ratio is! found to be: 
CReromec2 O30 To sbe=) 15414 
Thus the common divisor in the mole fraction calculation is 
given by: 
Chere Ge 6 /imet 30.6 oO) ete oa a2 oe eb AO) Un 8 au 
= ea PESTS 


iHeamoOle stractions are now foling from Eqs.) .4-D. 7 tosbe: 


Mee = 20.062) eatin Some o.UU 


MECN) PERSO CONT iss wen 4D > 


mf (CO) 5 23.0) Me OO t= eno, 
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Appendix E Computer Programs for Chapter 8 


This appendix contains the source code of the computer 
program which was used to integrate the six ordinary 
differential equation model presented in Chapter 8. Also 
included is the program which was used to determine the 
number, and stability, of steady-states for the six ODE 


model. 
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Cc 

C THIS PROGRAM IS USED TO INTEGRATE THE SYSTEM OF SIX ODE'S 
C PRESENTED IN CHAPTER 8. THE INTEGRATION CAN BE PERFORMED 
C USING BITHER GEAR'S METHOD OR THE RKF METHOD. 

Cc 

FTN4,X 


BLOCK DATA NAME, NAMED BLOCK DATA AREAS 

DOUBLE PRECISION F,K1,KM1,K2,KM2,K3,G1,G2,G3 
cee A2 AS, BT, B2,B3 

COMMON APARAMB F,K1,KM1,K2,KM2,K3,G1,G2,G3 
CAMA CA SeIB 1 Be BS 

COMMON /COUNT/ NPINT,IPINT,NDLAY,IMAX,IBUF(100) 
*, IDCBO( 144) ,VQO 

COMMON /MAN/ NINT,IPLOT,NAMEO(3) 

END 

PROGRAM DIM6(,31005), 6 ODE MODEL FOR CO OXIDATION 
EXTERNAL F6,PDERV 

DIMENSION X(7) 

DOUBLE PRECISION Y(6),TIME,H,HMX,AER,RER,DTIME 
*x,G(8,6) ,GSAVE(78) ,GMAX(6),GERR(6) ,GPW(72) 
COMMON /COUNT/ NPINT,IPINT,NDLAY,IMAX,I1BUF(100) 
* IDEBO (144) ;VO0 

COMMON /MAN/ NINT,IPLOT,NAMEO(3) 

CALLA DTACH 

CALL INPUT(LUNO,TIME,H,HMX,AER,RER,DTIME,N,Y) 
NPINT=1 

LF (NDDAY -GT.0) NPINT=O 

Om LO = 18 N 

GMAX(I )=0.01D0 

10 GA ia) = Vee) 

uS=0 

DORS0 SIPINT=a, NINT 

IF (INTG.EO. 0} 

*CALL RKF(TIME,Y,N,F6,DTIME,H,HMX,AER,RER,IFLAG) 
IF(INTG.EQ.1)CALL SGEAR(N,TIME,DTIME,G,GSAVE,H,HMX 
*,RER,1,GMAX,GERR,IFLAG,JS,6,GPW, PDERV) 

PeeCINTG. HOMO). GO T6730 

HOw 0 Li eiN 

20 VWORmEGGivr) 


¥ 


30 TECIFUAG INE S1))-GO TO-60 
PRINT eT INDUAY IGG TO. 56 
XT J=SNGU CPIME#VOO) 


DOWLOmLIe= 1 N 
40 X GDI+ PY=SNGLCECEE)) 
IF(IPINT/IPLOT*IPLOT.NE.IPINT) GO TO 50 
TAGE MA SENO. 1) GO T0250 
GALL? GODE 
WRITE (BUR VECO RAGA Ch 4% baer) 
800 FORMAT (12(1PE13.5)) 
CALL WRITF(IDCBO,IERR,IBUF, 46) 
WRITE (GUNG 4S 10 wes Cie Lae) 
810 FORMAT (7(1PE11.3)) 
NPINT=NPINT+1 
50 CONTINUE 
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WRITE(LUNO,820) NPINT 


820 BORMA TG" se OP TPOINTS =" 15) 
GOuTO WG 
60 WRITE(LUNO,830) IFLAG,TIME, (Y(I),1I=1,N) 
830 BORMA Tame TS 70D 15. 6.4 X)) 


70 CALL: ICHOSEX TDCBO) 

CALL OPEN(IDCBO,IERR,NAMEO,2,1,ICRO) 
CALL RWNDF(IDCBO) 
GALE CODE 
WRITE(IBUF,840) NPINT 

840 FORMAT (16) 
CALL WRITF(IDCBO,IERR,IBUF,3) 
GAEL CHOSE( IDCEO) 
STOP 
END 
SUBROUTINE INPUT(LUNO,TIME,H,HMX,AER,RER,DTIME,N,Y) 
DIMENSION X(7),NAMEI(3),IDCBI(144),IBUFA(50) ,IRBUF(33) 
DOUBLE PRECISION Y(6),TIME,H,HMX,AER,RER,DTIME 
DOUBLEMPRECLSTON “F Kt KMiIeK2.-KM2K3 Gi G2 "G37Aq Az A3 
Heap 2 Boe COn 07s COZ OS COS eT 
*,RK1,RK2,RK3,RKM1,RKM2,E1R,E2R,E3R,DH1,DH2,DH3 
eae Lm MASS. GP TO. POMOC, PCO.hO2, CONCT,DEXP,O 
COMMON /PARAM/ F,K1,KM1,K2,KM2,K3,G1,G2,G3 
Oa Ae A EB eB 2B 3 
COMMON /COUNT/ NPINT,IPINT,NDLAY,IMAX,IBUF(100) 
*, LDCBO( 144) VOO 
COMMON /MAN/ NINT,IPLOT,NAMEO(3) 
EQUIVALENCE (NAMEI(1),IRBUF(2)),(ICRI,IRBUF(6) ) 
EQUIVALENCE (ICRO,IRBUF(10)),(LUN,IRBUF(14)) 
LUNW=LOGLU(ISIS) 
LUNR=LUNW 
CALL IFMTE(LUNW) 
GAEL AGERSTCIBUFA,— 17 )lLOG) 
CALL PARSE(IBUFA,ILOG,IRBUF) 
PE CI HOG GT AGC Go, Tone 00 
WRITE (LUNW, 800) 

800 FORMAT(" NAME OF INPUT DATA FILE ") 
READ(LUNR,290) NAMEI 

290 FORMAT (3A2,1X,12) 
WRITE (LUNW, 801) 

801 FORMAT("ENTER INPUT CRN, OUTPUT CRN, AND OUTPUT LUN") 
READ(LUNR,*) ICRI,ICRO,LUN 

600 §LUNO=LUN 
CAL OPEN GI DCBIWI ERR NAMED 05.1, 1CRT) 
CALL READF(IDCBI,IERR,IBUF) 
CALIS CODE 
READ(IBUF,*) N,INTG,IDIMG,INTAL 
CALL READF(IDCBI,IERR,IBUF) 
GALE CODE 
READ(IBUF,*) TIME,DTIME,H,HMX,AER,RER 
CALL READF(IDCBI,IERR,IBUF) 
CALL CODE 
READ GU BURR eo (Oz, GO? AOS MCOS AT 
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CALL READF(IDCBI,IERR,IBUF) 

GAR EECODE 

READ(IBUF,*) RK1,RKM1,RK2,RKM2,RK3 
CALL READF(IDCBI,IERR,IBUF) 
GAMUGODE 

READ(IBUF,*) E1R,E2R,E3R,DH1,DH2,DH3 
CALL READF(IDCBI,IERR,IBUF) 
GAGTELGODE 

READ(IBUF,*) U,A,V,L,MASS,CP 
CALL READF(IDCBI,IERR,IBUF) 

CALI: «CODE 

READ(IBUF), *) TO,PO,Q0,FCO;FO2 
CONCT=P0/(8.314D0*TO) 

Q=Q0* (TO/273.D0)*(1.D+05/PO0) 
G1=E1R/TO 

G2=E2R/TO 

G3=E3R/TO 

K1=RK1*A/Q 

K2=RK2*A/Q 

K3=RK3*A*DEXP(-G3) /(Q*FCO*CONCT) 
KM1=RKM1*A*DEXP(-G1)/(Q*FCO*CONCT) 
KM2=RKM2*A*DEXP(-G2) /(Q*FO2*CONCT) 
B1=Q*DH 1*CONCT*FCO/(U*A*TO) 
B2=Q*DH2*CONCT*FO2/(UxA*TO) 
B3=Q*DH3 *CONCT*FCO/(U*XA*TO) 
Mies CONCTAYCO/ (ASL) 
A2=2.D0*V*CONCT#FO2/(A#L) 
A3=U*xA*V/(MASS*CPxQ) 

F=FCO+FO2 

VQO=V/0 

GO TO 30 ) 

CALL READF(IDCBI,IERR,IBUF) 

CALL CODE 

READ(IBUF,*) K1,KM1,K2,KM2,K3 
CALL READF(IDCBI,1IERR,IBUF) 

CAL, GODE 

READ(IBUF,*) G1,G2,G3,B1,B2,B3 
CRUG, READFUIDEBI® TRRR IBUP) 

CADIS CODE 

READ Gl BUF XIE AG, AZ, AS 
VOO=1.0D0 


CALL READF(IDCBI,IERR,IBUF) 


GARI CODE 

READ(IBUF,*) NINT,IPLOT,NDLAY 

PRGINTAT EO20)) ChEIIENTToCGOr O21, COS,0S) 
WRITE(LUNO,330) TIME,DTIME,H,HMX,AER,RER 
IF(IDIMG.EQ.0) GO TO 40 

WRITE(LUNO, 330) RK1,RK2,RK3,RKM1,RKM2 
WRITE(LUNO,330) E1R,E2R,E3R,DH1,DH2,DH3 
WROETE (LUNG? 330) 0U pAGV pb, MASS@CP 

WRITE (LUNO, 330) TO,PO,Q0,FCO,FO2 

WRITE GOUNG eis 0) tC OVO?) CO2FOS FCOS su 

WRCE CHUNG oS Uue fh TKM Is Ky KMee KS 

WRITE (LUNOW35I0) 5G 1, G27G3;Biaee Bs 
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WRETECLUNO}330)5R AT ,A2,A3 
eee FORMAT(12(1PD13.5)) 

WRITE(LUNO,345) NINT,IPLOT,N,INTG,NDLAY 
345 FORMAT (6(5X,14)) 


M@lr=CO 

V2 =O? 

Vis) =CO2 
V4 OS 
Vieow=COS 
Vicouvet 
X(1)=SNGL (TIME) 


BOM veri = 1,6 
1 re le SNGL Gy Gi),) 
CALL DFILE(IDCBI,ICRI,IDCBO,NAMEO,ICRO) 
CALL .CODE 
WRITE(IBUF,350)NINT 
350 FORMAT (16) 
CALL WRITF(IDCBO,IERR,IBUF,3) 
BeGBolNG SEO? 6) GOSTO 50 
CALL CODE 
WRITE(IBUF,360) RK1,RK2,RK3,RKM1,RKM2 
CALL WRITF(IDCBO,IERR,IBUF, 34) 
GAGE CODE 
WRITE(IBUF,360) E1R,E2R,E3R,DH1,DH2,DH3 
CALL WRITF(IDCBO,IERR,IBUF, 39) 
CARE VCODE 
WRITE(IBUF,360) U,A,V,L,MASS,CP 
CALL WRITF(IDCBO,IERR,IBUF, 39) 
CALE GODE 
WRITE(IBUF,360) TO,PO,Q0,FCO,FO2 
CALL WRITF(IDCBO,IERR,IBUF, 34) 
50 GALL SCODE 
WRITE(IBUF,360) K1,KM1,K2,KM2,K3 
360 FORMAT(12(1PE13.5)) 
CALL WRITF(IDCBO,IERR,IBUF, 34) 
CALE “CODE 
WRUTE CPBURY 260 )*G1-G27G37.Bi Be, BS 
CALL WRITF(IDCBO,IERR,IBUF,39) 
GAGE OC CODE 
WRITE( LBUF, S60)" F,A1,A2,A3 
CALL WRITF(IDCBO,IERR,IBUF, 26) 
| CALGSCODE 
WRITE(IBUF,360) TIME,DTIME,H,HMX,AER,RER 
CALL WRITF(IDCBO,IERR,IBUF, 39) 
Tr ONDEAY I GI.0))) SRETURN 
CADE ECORE 
WRITE GLBUF $360) CX(CT jeyta ie 
CALL WRITF(IDCBO,1ERR,IBUF,46) 
RETURN 
END 
SUBROUTENGSUENI TE (CO ,02,7T,COS, 0S) 
DOUBLE PRECISION. DEXP,DSORT,CO,COS,02,7T,OS,RM1, RICO, F 
#°RM2,R202,>K1)/KM1,K2,KM2,K3,G1,G2,G3,A1,A2,A3,81,32,B3 
COMMON /PARAM/ F,K1,KM1,K2,KM2,K3,G1,G2,G3 
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HPA, A2,AS Bly B2 B3 

RM1=KM1*DEXP(G1*T/(1.D0+T) ) 

RICOaK 1 «60 

RM2=KM2*DEXP(G2*T/(1.D0+T) ) 

RZO2Z=K2*O2 

66S=0..0D0 

OS=0.0D0 

Pea CO CET. 1.D=10)7AND,(O2. LT. 1.D-10) ) RETURN 
Nagar CO.. GH. 4 10) 2AND.(027GE ai. D- 10) )4G0: TO 20 
MeCEOVGEE i. D=10).GO: TO, 10 
OS=1.D0/(1.D0+DSORT(RM2/R202) ) 

RETURN 

COS="— DOV (4. DOFRMI/RICO) 

RETURN 
OS=1.D0/(1.D0+DSORT(RM2/R202)*(1.D0+R1CO/RM1) ) 
COS=1.D0-OS-DSORT(RM2/R202) *OS 

RETURN 

END 

SUBROUTINE F6(TIME,Y,YP) 

DOUBLE PRECISION Y(6),YP(6),TIME,DEXP 

DOUBLE MPRECISION ET, CO}02) CO2 (COS, 0S, RCO ,RCO2,7RO2, 000,0 
DOUBLE PREGISTONI DT , DEO, DO2,.DCO?2 ,DCOS , DOS, T1,S 
Pare Kehoe KM KS G1.G2,6G3.Al,A2,A3,B1,B2,B3 
COMMON /PARAM/ F,K1,KM1,K2,KM2,K3,G1,G2,G3 

CA ROT 3 Mp OBS 

COnY (a0) 

O2Z=Vi2) 

CO2=Y(3) 

OS=Y(4) 

COS= 705) 

fv 6) 

S=1.D0=-COS-OS 

T1=T+1.D0 

TecGun oO. 000) S=0. 000 
RCO=K1*CO*S-KM1*DEXP(G1*T/T1)*COS 
RO2=K2*02*S*S-KM2*DEXP(G2*T/T1) *OS*OS 
RCO2=K3*DEXP(G3#T/T1) *COS*0S 
Q00=1.D0+(F/(1.D0+2.D0*A1/A2) )*(2.D0*A1 

*x (RCO2-RCO) /A2-RO2) 

DCO=1.D0-QQ0*CO-RCO 

DO2=1.D0-000*02-RO2 


_ DCO2Z=-QQ0*CO2+RCO2 


DOS=A2*RO2-A1*RCO2 
DCOS=A1*(RCO-RCO2) 
DT=A3«(B1*RCO+B2*RO2+B3*RCO2-T) 
YP(A) =DCoO 

Me C2) =DOZ 

MES = DEOZ 

YP(4)=DOS 

YP.(5)=DCOS 

YP(6)=DT 

RETURN 

END 

SUBROUTINE PDERV(TIME,Y,YP,PW,M,IND) 
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DOUBLE PRECISION Y(8,M),YP(M),PW(M,M) , TIME,DEXP,DSQRT 


SDA Boy DA MOO, OF COZ COS 0S SO TIURT R2,R3-T 12 
Fa TeeDCOADOZsDCO2 DCOS, DOS; RCO, RCO? RO2 

Teh KM N2eKMey Ko Gt,G2;,G3,A1,A2,A3,B1,B2,B3 
COMMON /PARAM/ F,K1,KM1,K2,KM2,K3,G1,G2,G3 
APPA HAS AS BeBe, BS 

GO=VaGt 18) 


S=1.D0-COS-OS 

Teen ODO eS =O). 0D0 
Q0=00+(8.314D0#A*TO/P) *(-RCO-RO2+RCO2) 
RCO=K 1*CO*S-KM1*DEXP(G1*T/T1)*COS 
RO2=K2*02*S*S-KM2*DEXP(G2*T/T1) *OS*OS 
RCO2=K3*DEXP(G3*T/T1)*COS*0S 
GENGUND? BGeSW) <GOMTO: 100 
DCO=1.D0-F*CO-RCO 
DO2=1.D0-F*02-RO2 
DCO2=-F*CO2+RCO2 
DOS=A2*RO2-A1*RCO2 

DCOS=A1%* (RCO-RCO2) 
DT=A3%(B1*RCO+B2*RO2+B3*RCO2-T) 
VPC1) =DCO 

YP02)=DG2 

¥P(39)=DCO2 

YP(4)=DOS 

VP (> }=DCOS 

YP(6)=DT 

RETURN 

CONTINUE 

R1=KM1*DEXP(G1*T/T1) 
R2=KM2*DEXP(G2*T/T1) 
R3=K3*DEXP(G3*T/T1) 

ao eee 

PY GI = -RK ES 

PW 2) -0. 0D0 

pw(1,3)=0.0D0 


, PW(1,4)=K1*CO 
PW(1,5)=K1*CO+R1 
PW(1,6)=R1*G1*COS/T12 
PW(2,1)=0.0DO0 
PWU2, 2) =-P-K24S#S 
Pw(2,3)=0.0D0 
PW(2,4)=2.D0*R2*0S+2.D0*K2*O02*S 
PW(2,5)=2.D0*K2*02*S 
PW(2,6)=R2*G2*O0S*0S/T12 
Pw(3,1)=0.0D0 
PW(3,2)=0.0D0 
PW(3,3)=-F 
PW(3,4)=R3*COS 
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PW(3,5)=R3*0S 

PW(3,6)=R3%*G3*COS*0S/T12 

PW(4,1)=0.0D0 

PW 42) =A24K24*S4S 

PWi4) 3) =00D0 
PW(4,4)=-2.D0*A2*K2*02%*S-2.D0*A2*R2*OS-A1*R3%*COS 
PW(4,5)=-2.D0*A2*K2*02*S-A1*R3*0S 
PW(4,6)=-A2*R2#*G2*0S#0S/T12-A1#R3*G3*COS*0S/T12 
PW(52 J=A14KU eS 

PW(5,2)=0.0D0 

PWS 3) =00D0 

PW(5,4)=A1*(-K1*CO-R3*COS) 

ae 5)=A1*(-K1*CO-R3*OS-R1) 


W(5,6)=A1*(-R14G1#COS/T12-R3*G3%*0S*COS/T12) 

re Matar ics 

W(6,2)=A3*B2*K2*S*S 

PW(6,3)=0.0D0 

PW(6,4)=A3*(-B1*K1*CO-2.D0*B2% (K2*02*S+R2*0S ) 
*+B3*R3*COS) 

Wi 6, OJ =AG* (—B1#(K1*CO+R1)—2.D0#B2*K24*024S+B3*R3+*0S) 
PW(6,6)=A3% (-B1*R1*COS*G1/T12-B2#R2*G2*0S*0S/T12 
*+B3*R34*G3*0S*COS/T12-1.0D0) 

RETURN 

END 


This program also requires the RKF, SGEAR, MXINV and DFILE 
subroutines listed in Appendix B. 
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IMPLICIT REAL*8(A-H,K,0O-Z) 


PROGRAM UCALEJUATES @% MEY SrZ OS™ COS, T° GIVEN. ALD OTHER 
PARAMETERS FOR THE 6 ODE MODEL. AFTER CALCULATION OF 
THE ABOVE CONSTANTS, THE EIGENVALUES ARE CALULATED. 


EXTERNAL FCN 

DIMENSION EVAL(2,6) ,TOS(3,99) ,WK(68),TSS(99) ,OSS(99) 
mw COSS( 1OGN BRR (4100) -P(27) 

COMPLEX*16 W(6),ZSM,ZLG 


BOULVALENCE (WC i) ,BVAL Ci 1) 8(ZSM,0S7)_/(ZLG70S2) 
CP RR aGh (2) RKMI)), (PCS) RR2), (PCa) RRM2) ( PC5)7 
ERG) (PCUG)9 BIR), G27) E2R),(P(8) ,E3R),(P(9) ,DH1) 

PG LeDie) AUP él Lye Da) eb ute) | UL, GP CEs a) ACR IAD AD 
WAGES) eRE POE (Oly RMASS) (PGi) CP)  C(PC18 ) TO) 
ee iS ie DO ee ORK Oe OO 9 UP( 24) 7 FCO), (P0227 F02) 


COMMON T/P 0/2 ae 2, OS SCOS T 

COMMON /P2/ K1,KM1,K2,KM2,K3,A1,A2,A3 
GOMMON 9/237 G1,G2,/G3 ,B1,B2,B3;F 
GARMATRACER (at, 162)) 

READ(5,900) K1,KM1,K2,KM2,K3 

READ (5/900) FARPAZ)A3,B1,B2,B3 
READ(5;900) "G1,G2,G3,F 
FORMAT(9D13.5) 

READ(5,900) RK1,RKM1,RK2,RKM2,RK3 
READ(5,900) E1R,BE2R,E3R,DH1,DH2,DH3 
READ (5, 900} USA;V,RE, RMASS,CP 
see 900) TO, PO, OO -PCORPOCRE 
READ(5,905) NDATA 

FORMAT(14) 

DO 100 IDAT=1,NDATA 

READ(5,900) A3,B3 

READ(5,910) IPAR,PAR 
FORMAT(1I4,D13.6) 

P(IPAR)=PAR 

CONCT=P0/(8.314D0*TO) 

O=00* (T0/2732D0) 4 (1 .D+057P0) 
G1=E1R/TO 

G2=E2R/TO 

G3=E3R/TO 

K1=RK1*A/Q 
_K2=RK2*A/Q 

K3=RK3*A*DEXP(-G3) /(Q*FCO*CONCT ) 
KM 1=RKM1*A*DEXP(-G1) /(Q*FCO*CONCT ) 
KM2=RKM2*A*DEXP (-G2) /(Q*FO2*CONCT ) 
B1=Q*DH1*CONCT*FCO/(U*A*TO) 
B2=Q*DH2*CONCT*FO2/(U*A*TO) 
B3=Q*DH3*CONCT*FCO/(U*A*TO) 
A1=V*CONCT*FCO/(A*RL) 
A2=2.D0*V*CONCT#FO2/(A*RL) 
A3=U*xA*V/(RMASS*CP#Q) 

WRITE(6,805) 

FORMAT ( ann ad) ) 

WRITE(6, 806)RK1, RKM1,RK2,RKM2,RK3 
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FORMAT(10X,6(1PD11.4)) 

WRITE(6,806) E1R,E2R,E3R,DH1,DH2,DH3 
WRITE(6,806) U,A,V,RL,RMASS,CP 
WRITE(6,806) TO,PO,Q00,FCO,FO2 
WRITE(6,807) 


FORMAT(///) 

WRITE(6,810) K1,KM1,K3 

PORMA RCT PO Meer PD) 24° T30,K-1 = AIPDG2. 4 9T500N Kes 
Pe PD) 2; 43) 

WRITE(6,811) K2,KM2,F 

BORMAG UNTO; eK 2s =e PDO eT 30 8 K=2 =) opp 2 4 TSO. ver: 
wear Pl) 1 244.) 

WRIDEG6 RE 12) Al A2 AS 

BORMAWGnrOewsAveeee (PD OPA T3208 Age=' (PDP? 4oTS0 6YA Ss 
ie = Were oe.) 

WRETE (676 (3) 858 12 B253B3 

BORMAG (iO wea ts e CD (Onde Ts Oger eras Ppa. ANTS e Es | 
t= P24") 

WRUTE (ONG 2) 6C.1G2.G3 

BORMATmGTaO Gat 1 DDT Oars (Glas!) 1 PD ie 4650) “GS. 
eae PD 2s 4) 


TMAX=B1+B3+B2*A1/A2 

Done eb i= Ss 

POMC 2=4., 5 

PCSGie Cl teats onl) =i 1. ODOSIMAX 5. DO} +DFLOATCIZ—Ip 
+0.10*TMAX 

EOS ra) eS. EP) = 1 0D0/75. DOJ *£DFLOAT.CL t=1) 402 10D0 
POS Go) Gla) +5 nl2)=0.95D0-TOS (201 1-1) +5412) 

CONTINUE 

DORZ2Rio=t25 

NSIG=8 

N=3 

IMAX=75 

IER=0 

CALL UERSET(0,LEVOLD) 

HDOMCMeING=f33 

MOS CINGS Ig). =TOS CINC? 1a) 4100005 D0 . 
CALL ZSCNT(FCN,NSIG,N,IMAX,PAR,TOS(1,I1J),FNORM,WK,IER) 
DOMZGHINC=1 "3 

TOSGINC PO y=TOS GING Fla) mI 000, D0 

IERR(IJ)=IER 


_NSS=0 


BO) Puan ame ees 

TeCHMRR Ol BO vic S260 shOez2 So 

Te (@TOS GIS EI NEE f02) OR TOS (1 1d) GE iMAR eGo ntC mes 
GE GGTOS CONS) ZERTOT) ORS CTOS (2 (LJ) GE SoD OUe homme. > 
IF(NSS.EQ.0) GO TO 24 

WOmls PK = NSS 

REL=1.D-06 

tet (DABS (@CTOS (1, 10)-TSS (TK) )7TSS(IK IGT REL) SANDE 
G(DABS (TOS 62°] ) —COSS (1K) )/COSS CLK) ) SLT. REL) eGo. FOP 25 
CONTINUE 

NSS=NSS+1 

™SS (NSS }=TOS(4, 10) 
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€OSS(NSS)=TOS (2970) 
SSSuNSS ieTOS03 1 J) 
GOLTOR25 
24 CONTINUE 
NSS=1 
BSSUNSS p=TOSW lu 
COSS NSS] ROS (210) 
GSSENSS)=TOS (3. Td) 
25 CONTINUE 
DO Z0SI2=fSnss 
MeTSS Ord) 
GOS=coss(1 J) 
OS=0SS(2J} 
%=(1,D0-T/TMAX) /P 
¥=(1.D0-A1i*¥T/ (AZ2*TMAX))/F 
Z=K3*DEXP(G3*T/(1.D0+T) )*OS*COS/F 
FP(COSTEOS SCR 1201D0).0R: (OS+COS (ET -0201D0)3GO To 30 
PaCS + 2G ley eOeDo) OR, kt 2 Um. 0. S8D0I)— GO TO 30 
CALL EVALUE(W) 
WRITE COMB 260" 1, % iy 


830 FORMAT (TS, 6 

WRITE(6, 830 
30 CONTINUE 
100 CONTINUE 
200 CONTINUE 

STOP 

END 

SUBROUTINE EVALUE(W) 

IMPLICIT REAL*¥8(A-H,K,O-Z) 

DIMENSION RJ(6,6),WK(6) ,BB(6) 

COMPLEX*16 W(6),2Z2Z(6,6) 

COMMOND PU REYeZF0S, COSer 

COMMON /P2/ K1,KM1,K2,KM2,K3,A1,A2,A3 

COMMON /P3/ G1,G2,G3,B1,B2,B3,F 

IER=0 

=i 2D0-COS-OS 
“Ttets 1 7 D0 

aoe | 

R1=KM1*DEXP(G1*T/T1) 

R2=KM2*DEXP(G2*T/T1) 

R3=K3*DEXP(G3*T/T1) 


826 ea ae SMe PD. 6, 130, Kee PD oo. 
<a> Cee =P We 66) 
Beanies 827) Z5OS, COS 
827 FORMAT(T10, ve bee PD 6 o 0) OS= 8) 1PDis 6 roo 
eCOS= AEDs 67) 
WRITE(6, SSC MECEVAD C17 J)), i=, 69 
Ci 
) 


( 
1 POWs Rex iy 
(CN ie (20rd lees =a Gn) 


RIC LySSESK1eS 
Roi 2) OC oe 0 

Rot 33) = 000 

Ro (1, S)=Kiex 

RIC, S)=KhTeetR1 
RJ(1,6)=RU#G t*COS/T12 
Rapa.) = Cre 0 


a % Lg 
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RJ(2,2)=-F-K2*S*S 

Re 2 2 }=O08-0D0 
RJ(2,4)=2.D0#R2*OS+2.D0*K2*YxS 

Ree a 5= 2, DOeK2Z4V ES 

RU-G2 pow R24G2ZFOS*OS 7 T 12 

Ro €2 7 b)y=070D0 

Rowse) =0 20D0 

Ra (3) 0) =-F 

Ras 4). =R3*COS 

Rute >) =RSECS 

Renter p= RaxG 3 4COS*OS/T12 

Ravaei) =O 0D0 

Ro (42) =A2eK2*S4S 

RJ (4,3)=0.0D0 

RJ(4,4)=-2.D0*A2* (K2*Y*S+R2*OS )-A1*R3*COS 

RJ (4,5)=-2.D0*#A2*K2*VY*S-A1#R3*0S 
RJ(4,6)=-A2*R2*G2#0S*0S/T12-A1*R3*G3*COS*O0S/T12 
Re 5) th) =A eK es, 

Rd Gore = 0 0D0 

Ro5 7.3) =0. 0D0 

RU (5,4) =-A1*(K1*K+R3*COS) 
RJ(5,5)=-A1*(K1*X+R1+R3*0S) 

ROG .G l=-Ads CR 1 4G) *6GOS/T 12+R3*634COS*0S/T 12) 
RUA 6, 1) =AStB1eK 1*S 

RoCG 2) =ASEER ASK 24S eS 

Re 66) =020D0 

RJ(6,4)=A3*(-B1*K1*X-2.D0#B2% (K2*Y*S+R2*OS) +B3*R3*COS) 


RJ(6,5)=A3*(-B1*(K1#*X+R1)-2.D0*B2*K2*Y*S+B3*R3%0S) 


RJ(6,6)=A3*(-B1*R1*G1*COS/T12-B2*R2*G2*0S*0S/T12 


*+B3*R3*G3*O0S*COS/T12-1.D0) 


GALE EBGRE(RI,6,6,0,W, 22,6, WK, 1ER) 

TP CGLER TET. 32) RETURN 

WRITE(6,810) IER 

FORMAT(//,T20,'IER FROM EIGRF = ',14) 
RETURN 

END 

SUBROUTINE FCN(XX,FF,N, PAR) 

IMPLICIT REAL*8(A-H,K,O-Z) 

DIMENSION XX(N),FF(N),PAR(1) 

COMMON /P2/ K1,KM1,K2,KM2,K3,A1,A2,A3 
COMMON /P3/ G1,G2,G3,B1,B2,B3,F 
TMAX=B1+B3+B2*A1/A2 


ie XXL hoe CDC ek ay =O DU 


De OCC) (GRA TMAX# 100002 DO) XX (1 )J=TMAX 41 OOU Gee 
Pex) (0000. BO 
COS=ckt2 e000 D0 
OS=xKX(3)/10000.D0 
F{=DEXP(G1*#T/(1.D0+T) 
F2=DEXP(G2*T/(1.D0+T) 
X=(1.DO0+KM1*F1*COS) /( 
Y=(1.DO0+KM2*F2*x0S*0OS) 
ei) Glee DO) 
S=1.D0-COS-OS 
F3=DEXP(G3T1) 


) 

) 

F+K1*(1.D0-OS-COS) ) 

/ (B+K2%(-C14 D0=0S-GOS ) + x2oe) 
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FF(1)=B1*(1.D0-FP*X)+B2*(1.D0-F*Y)+B3*K3*F3*COS*O0S-T 
FF(2)=K1*X*S-KM1*DEXP(G1#T1) *COS-K3*DEXP(G3*T1) *COS#0S 
FF (3)=A2% (K2*Y*S*S-KM2*DEXP(G2*T1)*OS*OS)-A1*K3*F3 
**xCOS*OS 

Bet) =PPi7)410000-D0 

FE(2)=FP(2)*100,D0 

FE (3)=FF(3)4#100.D0 

RETURN 

END 
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